Coverage for gamdpy/calculators/calculator_hydrodynamic_correlations.py: 73%

1import numpy as np 

2import gamdpy as gp 

3from numba import jit 

4import math 

5import cmath 

6 

7@jit(nopython=True) 

8def __store__(dk, dmk, jk, jmk, sampletype, mass, pos, vel, ptypes, npart, index, nwaves, lbox): 

9 

10 I = complex(0.0, 1.0) 

11 

12 for k in range(nwaves): 

13 dk[index, k], dmk[index, k] = complex(0.0, 0.0), complex(0.0, 0.0) 

14 jk[index, k], jmk[index, k] = complex(0.0, 0.0), complex(0.0, 0.0) 

15 

16 kwave = 2.0*math.pi*(k+1)/lbox 

17 

18 for n in range(npart): 

19 if ptypes[n] == sampletype or sampletype == -1: 

20 kfac = cmath.exp(I*kwave*pos[n]) 

21 mkfac = cmath.exp(-I*kwave*pos[n]) 

22 

23 dk[index, k] = dk[index, k] + mass[n]*kfac 

24 dmk[index, k] = dmk[index, k] + mass[n]*mkfac 

25 

26 jk[index, k] = jk[index, k] + mass[n]*vel[n]*kfac 

27 jmk[index, k] = jmk[index, k] + mass[n]*vel[n]*mkfac 

28 

29 

30@jit(nopython=True) 

31def __calccorr__(dacf, jacf, dk, dmk, jk, jmk, nwaves, lvec): 

32 

33 for k in range(nwaves): 

34 for n in range(lvec): 

35 for nn in range(lvec-n): 

36 dacf[n, k] = dacf[n,k] + dk[nn,k]*dmk[nn+n, k] 

37 jacf[n, k] = jacf[n,k] + jk[nn,k]*jmk[nn+n, k] 

38 

39 

40 

41class CalculatorHydrodynamicCorrelations: 

42 """ 

43 Calculates the tranverse current autocorrelation function (jacf) and  

44 the longitudinal density autocorrelation function (dacf).  

45 

46 Example: See hydrocorr.py 

47 

48 Initialization variables 

49 - configuration: Instance of the configuration class 

50 - dtsample: Time between sampling (= integrator time step * steps_per_timeblock) 

51 - ptype: The particle type for which the correlations are calculated (default 0)  

52 - nwaves: Number of wavevectors (default 10) 

53 - lvec: Length of the time vector - sample time span then dtsample*lvec (defulat 100)  

54 - verbose: Write a bit to screen (default False) 

55 

56 Output:  

57 With the method read() data is written to files and returned to user  

58 """ 

59 

60 def __init__(self, configuration, dtsample, ptype=0, nwaves=10, lvec=100, verbose=False): 

61 self.conf = configuration 

62 self.ptype = ptype 

63 self.nwaves = nwaves 

64 self.lvec = lvec 

65 self.dt = dtsample 

66 

67 self.nsample = 0 

68 self.index = 0 

69 

70 self.lbox = configuration.simbox.get_lengths()[1] 

71 

72 # Storage arrays 

73 self.dk = np.zeros( (lvec, nwaves), dtype=np.complex64) 

74 self.dmk = np.zeros( (lvec, nwaves), dtype=np.complex64) 

75 

76 self.jk = np.zeros( (lvec, nwaves), dtype=np.complex64) 

77 self.jmk = np.zeros( (lvec, nwaves), dtype=np.complex64) 

78 

79 # Correlation array - density (longitudinal) 

80 self.dacf = np.zeros( (lvec, nwaves), dtype=np.complex64) 

81 

82 # Current density acf (transverse) 

83 self.jacf = np.zeros( (lvec, nwaves), dtype=np.complex64) 

84 

85 

86 if verbose: 

87 print(f"Types {self.ptype}, no. wavevectors {self.nwaves}, vector length {self.lvec}") 

88 

89 

90 

91 

92 def update(self): 

93 

94 __store__(self.dk, self.dmk, self.jk, self.jmk, self.ptype, 

95 self.conf['m'][:], self.conf['r'][:,0], self.conf['v'][:,1], 

96 self.conf.ptype[:], self.conf.N, self.index, self.nwaves, self.lbox) 

97 

98 self.index = self.index + 1 

99 

100 if self.index == self.lvec: 

101 __calccorr__(self.dacf, self.jacf, self.dk, self.dmk, self.jk, self.jmk, self.nwaves, self.lvec) 

102 

103 self.nsample = self.nsample + 1 

104 self.index = 0 

105 

106 

107 

108 

109 def read(self, save=True, fname_dacf="dacf.dat", fname_jacf="jacf.dat"): 

110 

111 if self.nsample > 0: 

112 

113 volume = self.conf.simbox.get_volume() 

114 

115 

116 out_dacf = np.zeros( (self.lvec, self.nwaves) ) 

117 out_jacf = np.zeros( (self.lvec, self.nwaves) ) 

118 

119 file_dacf = open(fname_dacf, "w") 

120 file_jacf = open(fname_jacf, "w") 

121 

122 for n in range(self.lvec): 

123 

124 fac = 1.0/(self.nsample*(self.lvec-n)*volume) 

125 

126 file_dacf.write("%f " % (self.dt*n)) 

127 file_jacf.write("%f " % (self.dt*n)) 

128 

129 for k in range(self.nwaves): 

130 out_dacf[n,k] = self.dacf[n, k].real*fac 

131 out_jacf[n,k] = self.jacf[n, k].real*fac 

132 file_dacf.write("%f " % (out_dacf[n,k])) 

133 file_jacf.write("%f " % (out_jacf[n,k])) 

134 

135 file_dacf.write("\n") 

136 file_jacf.write("\n") 

137 

138 file_dacf.close() 

139 file_jacf.close() 

140 

141 return (out_dacf, out_jacf) 

142 

143 else: 

144 return (0,0)