1module TMDGenData !********************************************************************************* 2! 3! Common data for TMDGen 4! 5!--------------------------------------------------------------------------------------------------- 6! 7! Intel Fortran 8! 9! Alexey N. Volkov, University of Alabama, avolkov1@ua.edu, Version 13.00, 2020 10! 11!--------------------------------------------------------------------------------------------------- 12 13use TPMGeom 14 15implicit none 16 17!--------------------------------------------------------------------------------------------------- 18! Constants 19!--------------------------------------------------------------------------------------------------- 20 21 integer*4, parameter :: MAX_TUBE = 1000000 ! Maximum number of tubes in 3D 22 23 real*8, parameter :: K_MDDU = K_MDMU / K_MDLU / K_MDLU / K_MDLU ! MD density unit (kg/m**3) 24 25 ! 26 ! These parameters are specific for carbon nanotubes and taken from module TubePotBase 27 ! 28 29 real*8, parameter :: TPbConstD = 5.196152422706632d+00 ! = 3.0**1.5 30 ! Mass of C atom 31 real*8, parameter :: TPBM = 12.0107d+00 ! (a.m.u.) 32 ! Lattice parameter and numerical density of atoms for a graphene sheet, see Dresselhaus et al, Carbon 33(7), 1995 33 real*8, parameter :: TPBA = 1.421d+00 ! (Angstrom) 34 real*8, parameter :: TPBD = 4.0d+00 / ( TPBConstD * TPBA * TPBA ) ! (1/Angstrom^2) 35 ! Specific heat of carbon nanotubes 36 real*8, parameter :: TPBSH = 600.0d+00 / K_MDCU ! (eV/(Da*K)) 37 38!--------------------------------------------------------------------------------------------------- 39! Governing parameters 40!--------------------------------------------------------------------------------------------------- 41 42 ! Parameters of the sample 43 44 real*8 :: LS0 = 4000.0 ! Sample size in x, y-directions (Angstrom) 45 real*8 :: HS0 = 4000.0 ! Sample size in z-direction (Angstrom) 46 real*8 :: DS0 = 0.01 ! Density (g/cm**3) 47 integer*4 :: BC_X0 = 1 ! Boundary conditions in x-direction: 0, free; 1, periodic 48 integer*4 :: BC_Y0 = 1 ! Boundary conditions in y-direction: 0, free; 1, periodic 49 integer*4 :: BC_Z0 = 1 ! Boundary conditions in z-direction: 0, free; 1, periodic 50 51 ! Parameters of tubes 52 integer*4 :: ChiIndM = 10 ! Chirality index m of nanotubes 53 integer*4 :: ChiIndN = 10 ! Chirality index n of nanotubes 54 real*8 :: LT0 = 2000.0 ! Characterstic length of tubes (Angstrom) 55 integer*4 :: SegType = 0 ! 0, number of segments per tube is fixed 56 ! 1, rounded length of segments is fixed 57 integer*4 :: NSeg0 = 100 ! Number of segments per tube 58 real*8 :: LSeg0 = 20.0d+00 ! Length of the segment (Angstrom) 59 60 61 ! Parameters controlling the sample structure 62 63 real*8 :: DeltaT = 3.0 ! Minimal distance between tubes (Angstrom) 64 integer*4 :: NAmax = 50000 ! Maximal number of attempts (for SampleType = 4 it is used as an input paramtere for number of tubes) 65 real*8 :: GeomPrec = 1.0d-06 ! Geometrical precision 66 67!--------------------------------------------------------------------------------------------------- 68! Computed data 69!--------------------------------------------------------------------------------------------------- 70 71 real*8 :: RT0 = 6.785 ! Radius of tubes (Angstrom) 72 73 real*8 :: VS0 ! Desired volume of the sample, Angstrom**3 74 real*8 :: MS0 ! Desired mass of the sample, Da (For SampleType = 4 it is the defined fixed mass- definition is given in TMDGen7T) 75 76 real*8 :: CTCD ! Center to center distance between any surrounding tube and center tube (used for SampleType == 4 only) 77 78 integer*4 :: NT ! Real number of tubes 79 real*8, dimension(0:MAX_TUBE-1) :: RT ! Radii of tubes, Angstrom 80 real*8, dimension(0:MAX_TUBE-1) :: LT ! Lengths of tubes, Angstrom 81 real*8, dimension(0:MAX_TUBE-1,0:2) :: CT ! Coordinates of tubes' centers, Angstrom 82 real*8, dimension(0:MAX_TUBE-1,0:2) :: DT ! Directions of tubes 83 integer*4, dimension(0:MAX_TUBE-1) :: AT ! Parent axes of tubes. It is used only in GeneratorBundle () 84 85contains !****************************************************************************************** 86 87!--------------------------------------------------------------------------------------------------- 88! Pseudo-random number generator 89!--------------------------------------------------------------------------------------------------- 90 91 real*8 function randnumber () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 92 ! This function returns a pseudo-random number with uniform distribution in [0,1] 93 !------------------------------------------------------------------------------------------- 94 call random_number ( randnumber ) 95 end function randnumber !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 96 97 subroutine SetRandomSeed () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 98 ! This subroutine sets random seed for the pseudo-random number generator 99 !------------------------------------------------------------------------------------------- 100 integer :: i, n, clock 101 integer, dimension(:), allocatable :: seed 102 !------------------------------------------------------------------------------------------- 103 call RANDOM_SEED ( size = n ) 104 allocate ( seed(n) ) 105 call SYSTEM_CLOCK ( COUNT = clock ) 106 seed = clock + 37 * (/ (i - 1, i = 1, n) /) 107 call RANDOM_SEED ( PUT = seed ) 108 deallocate ( seed ) 109 end subroutine SetRandomSeed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 110 111!--------------------------------------------------------------------------------------------------- 112! Generators for (random) properties of nanotubes 113!--------------------------------------------------------------------------------------------------- 114 115 real*8 function TubeMass ( i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 116 ! This function returns the mass of the tube in Da 117 !------------------------------------------------------------------------------------------- 118 integer*4, intent(in) :: i 119 !------------------------------------------------------------------------------------------- 120 TubeMass = M_2PI * RT(i) * LT(i) * TPBM * TPBD 121 end function TubeMass !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 122 123 real*8 function TubeSpecificHeat ( i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 124 ! This function returns the specific heat of the tube 125 !------------------------------------------------------------------------------------------- 126 integer*4, intent(in) :: i 127 !------------------------------------------------------------------------------------------- 128 TubeSpecificHeat = TPBSH 129 end function TubeSpecificHeat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 130 131!--------------------------------------------------------------------------------------------------- 132! Reading and printing of input parameters 133!--------------------------------------------------------------------------------------------------- 134 135 subroutine LoadGoverningParameters () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 136 ! This function reads governing parameters from xdt file 137 !------------------------------------------------------------------------------------------- 138 integer*4 :: Fuid, i 139 character*512 :: Msg 140 !------------------------------------------------------------------------------------------- 141 Fuid = OpenFile ( 'TMDGen.xdt', 'rt', '' ) 142 read ( unit = Fuid, fmt = '(e22.12)' ) LS0 143 read ( unit = Fuid, fmt = '(e22.12)' ) HS0 144 read ( unit = Fuid, fmt = '(e22.12)' ) DS0 145 read ( unit = Fuid, fmt = '(i22)' ) BC_X0 146 read ( unit = Fuid, fmt = '(i22)' ) BC_Y0 147 read ( unit = Fuid, fmt = '(i22)' ) BC_Z0 148 read ( unit = Fuid, fmt = '(i22)' ) ChiIndM 149 read ( unit = Fuid, fmt = '(i22)' ) ChiIndN 150 read ( unit = Fuid, fmt = '(e22.12)' ) LT0 151 read ( unit = Fuid, fmt = '(i22)' ) SegType 152 read ( unit = Fuid, fmt = '(i22)' ) NSeg0 153 read ( unit = Fuid, fmt = '(e22.12)' ) LSeg0 154 read ( unit = Fuid, fmt = '(e22.12)' ) DeltaT 155 read ( unit = Fuid, fmt = '(i22)' ) NAmax 156 read ( unit = Fuid, fmt = '(e22.12)' ) GeomPrec 157 call CloseFile ( Fuid ) 158 end subroutine LoadGoverningParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 159 160 subroutine PrintGoverningParameters () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 161 ! This function prints governing parameters to xlg file 162 !------------------------------------------------------------------------------------------- 163 integer*4 :: Fuid, i 164 !------------------------------------------------------------------------------------------- 165 Fuid = OpenFile ( 'TMDGen.xlg', 'wt', '' ) 166 write ( unit = Fuid, fmt = '(e22.12,a)' ) LS0, ' : LS0, Angstrom' 167 write ( unit = Fuid, fmt = '(e22.12,a)' ) HS0, ' : HS0, Angstrom' 168 write ( unit = Fuid, fmt = '(e22.12,a)' ) DS0, ' : DS0, g/cm**3' 169 write ( unit = Fuid, fmt = '(e22.12,a)' ) DS0, ' : SC0, 1/A**2' 170 write ( unit = Fuid, fmt = '(i22,a)' ) BC_X0, ' : BC_X0' 171 write ( unit = Fuid, fmt = '(i22,a)' ) BC_Y0, ' : BC_Y0' 172 write ( unit = Fuid, fmt = '(i22,a)' ) BC_Z0, ' : BC_Z0' 173 write ( unit = Fuid, fmt = '(i22,a)' ) ChiIndM, ' : ChiIndM' 174 write ( unit = Fuid, fmt = '(i22,a)' ) ChiIndN, ' : ChiIndN' 175 write ( unit = Fuid, fmt = '(e22.12,a)' ) LT0, ' : LT0, Angstrom' 176 write ( unit = Fuid, fmt = '(i22,a)' ) SegType, ' : SegType' 177 write ( unit = Fuid, fmt = '(i22,a)' ) NSeg0, ' : NSeg0' 178 write ( unit = Fuid, fmt = '(e22.12,a)' ) LSeg0, ' : LSeg0, Angstrom' 179 write ( unit = Fuid, fmt = '(e22.12,a)' ) DeltaT, ' : DeltaT' 180 write ( unit = Fuid, fmt = '(i22,a)' ) NAmax, ' : NAmax' 181 write ( unit = Fuid, fmt = '(e22.12,a)' ) GeomPrec, ' : GeomPrec' 182 call CloseFile ( Fuid ) 183 end subroutine PrintGoverningParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 184 185!--------------------------------------------------------------------------------------------------- 186! Printing of sample parameters 187!--------------------------------------------------------------------------------------------------- 188 189 subroutine PrintSampleParameters ( ParType ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 190 ! This function prints the most imprtant parameters of the sample. 191 ! In the code, it used twice to print parameters of the desired and really generated samples. 192 !------------------------------------------------------------------------------------------- 193 character*(*), intent(in) :: ParType 194 real*8 :: MP, M, V 195 !------------------------------------------------------------------------------------------- 196 print '(a,a,a)', '*** ', trim(ParType), ' properties of the sample' 197 print '(a34,a,f15.4,a)', 'L', ' : ', LS0, ' A' 198 print '(a34,a,f15.4,a)', 'H', ' : ', HS0, ' A' 199 print '(a34,a,f15.4,a)', 'Density', ' : ', DS0, ' g/cm**3' 200 print '(a34,a,e15.8,a)', 'Volume', ' : ', VS0, ' A*3' 201 print '(a34,a,e15.8,a)', 'Mass', ' : ', MS0, ' Da' 202 print '(a34,a,i10)', 'BC_X', ' : ', BC_X0 203 print '(a34,a,i10)', 'BC_Y', ' : ', BC_Y0 204 print '(a34,a,i10)', 'BC_Z', ' : ', BC_Z0 205 end subroutine PrintSampleParameters !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 206 207!--------------------------------------------------------------------------------------------------- 208! Initializing of basic geometrical parameters of the generated sample 209!--------------------------------------------------------------------------------------------------- 210 211 subroutine InitSample () !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 212 ! This function initializes the geometrical parameters of the sample (sizes, etc.) 213 !------------------------------------------------------------------------------------------- 214 215 BC_X = BC_X0 216 BC_Y = BC_Y0 217 BC_Z = BC_Z0 218 DomXmin = - LS0 / 2.0d+00 219 DomXmax = LS0 / 2.0d+00 220 DomYmin = - LS0 / 2.0d+00 221 DomYmax = LS0 / 2.0d+00 222 DomZmin = - HS0 / 2.0d+00 223 DomZmax = HS0 / 2.0d+00 224 225 if ( BC_X0 == 0 ) then 226 DomXmin = 0.0d+00 227 DomXmax = LS0 228 end if 229 if ( BC_Y0 == 0 ) then 230 DomYmin = 0.0d+00 231 DomYmax = LS0 232 end if 233 if ( BC_Z0 == 0 ) then 234 DomZmin = 0.0d+00 235 DomZmax = HS0 236 end if 237 238 DomLX = DomXmax - DomXmin 239 DomLY = DomYmax - DomYmin 240 DomLZ = DomZmax - DomZmin 241 DomLXHalf = 0.5d+00 * DomLX 242 DomLYHalf = 0.5d+00 * DomLY 243 DomLZHalf = 0.5d+00 * DomLZ 244 245 DS0 = DS0 / ( K_MDDU / 1.0d+03 ) 246 VS0 = LS0 * LS0 * HS0 247 MS0 = DS0 * VS0 248 end subroutine InitSample !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 249 250!--------------------------------------------------------------------------------------------------- 251! A few auxiliary functions 252!--------------------------------------------------------------------------------------------------- 253 254 subroutine GetTubeEnds ( X0, X1, i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 255 ! This function calculates coordinates of two ends of nanotube i 256 !------------------------------------------------------------------------------------------- 257 real*8, dimension(0:2), intent(out) :: X0, X1 258 integer*4, intent(in) :: i 259 !------------------------------------------------------------------------------------------- 260 real*8 :: LT2 261 !------------------------------------------------------------------------------------------- 262 LT2 = 0.5d+00 * LT(i) 263 X0 = CT(i,0:2) - LT2 * DT(i,0:2) 264 X1 = CT(i,0:2) + LT2 * DT(i,0:2) 265 end subroutine GetTubeEnds !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 266 267 logical function IsTubeInside ( i ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 268 ! This function returns true if nanotube i lies inside the sample. Otherwise it returns false. 269 !------------------------------------------------------------------------------------------- 270 integer*4, intent(in) :: i 271 !------------------------------------------------------------------------------------------- 272 integer*4 :: n 273 real*8, dimension(0:2) :: X0, X1, Xmin, Xmax 274 !------------------------------------------------------------------------------------------- 275 IsTubeInside = .true. 276 if ( BC_X == 1 .and. BC_Y == 1 .and. BC_Z == 1 ) return 277 call GetTubeEnds ( X0, X1, i ) 278 do n = 0, 2 279 Xmin(n) = min ( X0(n), X1(n) ) 280 Xmax(n) = max ( X0(n), X1(n) ) 281 end do 282 IsTubeInside = .false. 283 if ( BC_X == 0 .and. ( Xmin(0) < DomXmin .or. Xmax(0) > DomXmax ) ) return 284 if ( BC_Y == 0 .and. ( Xmin(1) < DomYmin .or. Xmax(1) > DomYmax ) ) return 285 if ( BC_Z == 0 .and. ( Xmin(2) < DomZmin .or. Xmax(2) > DomZmax ) ) return 286 IsTubeInside = .true. 287 end function IsTubeInside !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 288 289end module TMDGenData !***************************************************************************** 290