TABLE OF CONTENTS
- ABINIT/m_pimd
- m_pimd/pimd_apply_constraint
- m_pimd/pimd_coord_transform
- m_pimd/pimd_destroy
- m_pimd/pimd_diff_stress
- m_pimd/pimd_energies
- m_pimd/pimd_force_transform
- m_pimd/pimd_forces
- m_pimd/pimd_init
- m_pimd/pimd_init_qtb
- m_pimd/pimd_initvel
- m_pimd/pimd_is_restart
- m_pimd/pimd_langevin_forces
- m_pimd/pimd_langevin_random
- m_pimd/pimd_langevin_random_bar
- m_pimd/pimd_langevin_random_init
- m_pimd/pimd_langevin_random_qtb
- m_pimd/pimd_mass_spring
- m_pimd/pimd_noseehoover_forces
- m_pimd/pimd_nosehoover_propagate
- m_pimd/pimd_nullify
- m_pimd/pimd_predict_taylor
- m_pimd/pimd_predict_verlet
- m_pimd/pimd_print
- m_pimd/pimd_skip_qtb
- m_pimd/pimd_stresses
- m_pimd/pimd_temperature
- m_pimd/pimd_type
ABINIT/m_pimd [ Modules ]
NAME
m_pimd
FUNCTION
This module provides several routines and datatypes for the Path-Integral Molecular Dynamics (PIMD) implementation.
COPYRIGHT
Copyright (C) 2010-2024 ABINIT group (GG,MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_pimd 24 25 use defs_basis 26 use m_abicore 27 use m_dtset 28 use m_errors 29 use m_io_tools 30 use m_random_zbq 31 32 use m_numeric_tools, only : uniformrandom 33 use m_symtk, only : matr3inv 34 use m_geometry, only : mkradim 35 36 implicit none 37 38 private 39 40 !public procedures 41 public :: pimd_init 42 public :: pimd_nullify 43 public :: pimd_destroy 44 public :: pimd_init_qtb 45 public :: pimd_skip_qtb 46 public :: pimd_print 47 public :: pimd_is_restart 48 public :: pimd_temperature 49 public :: pimd_initvel 50 public :: pimd_langevin_random 51 public :: pimd_langevin_random_qtb 52 public :: pimd_langevin_random_bar 53 public :: pimd_langevin_random_init 54 public :: pimd_energies 55 public :: pimd_forces 56 public :: pimd_langevin_forces 57 public :: pimd_nosehoover_forces 58 public :: pimd_stresses 59 public :: pimd_diff_stress 60 public :: pimd_predict_taylor 61 public :: pimd_predict_verlet 62 public :: pimd_nosehoover_propagate 63 public :: pimd_coord_transform 64 public :: pimd_force_transform 65 public :: pimd_apply_constraint 66 public :: pimd_mass_spring
m_pimd/pimd_apply_constraint [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_apply_constraint
FUNCTION
Modify forces to take into account an holonomic constraint according to "pimd_constraint" parameter Available constraints: 0: no constraint 1: linear combination of coordinates
INPUTS
constraint=type of constraint to be applied mass(natom)=fictitious masses of atoms natom=number of atoms trotter=Trotter number wtatcon(3,natom)=weights for atomic constraints xcart(3,natom,trotter)=cartesian coordinates of atoms
OUTPUT
constraint_output(2)=several (real) data to be output when a constraint has been applied
SIDE EFFECTS
forces(3,natom,trotter)=array containing forces
SOURCE
2367 subroutine pimd_apply_constraint(constraint,constraint_output,forces,mass,natom,& 2368 & trotter,wtatcon,xcart) 2369 2370 !Arguments ------------------------------------ 2371 !scalars 2372 integer,intent(in) :: constraint,natom,trotter 2373 !arrays 2374 real(dp),intent(in) :: mass(natom),wtatcon(3,natom),xcart(3,natom,trotter) 2375 real(dp),intent(out) :: constraint_output(2) 2376 real(dp),intent(inout) :: forces(3,natom,trotter) 2377 !Local variables------------------------------- 2378 !scalars 2379 integer :: iatom,ii,iimage 2380 real(dp) :: af,lambda_cst,masstot,one_over_trotter,xcart_centroid,zz 2381 !character(len=500) :: msg 2382 !arrays 2383 real(dp) :: force_centroid(3),lambda_com(3),mat(3,3),matinv(3,3),vec(3),weightsum(3) 2384 2385 !************************************************************************ 2386 2387 2388 select case(constraint) 2389 2390 !=== No constraint ======================================================= 2391 case(0) 2392 2393 constraint_output(:)=zero 2394 return 2395 2396 !=== Linear combination of centroid coordinates ========================== 2397 case(1) 2398 2399 ! Some useful quantities 2400 masstot=sum(mass)*dble(trotter) 2401 weightsum(:)=sum(wtatcon,dim=2)*dble(trotter) 2402 zz=zero;af=zero;force_centroid=zero 2403 do iimage=1,trotter 2404 do iatom=1,natom 2405 do ii=1,3 2406 force_centroid(ii)=force_centroid(ii)+forces(ii,iatom,iimage) 2407 af=af+wtatcon(ii,iatom)*forces(ii,iatom,iimage)/mass(iatom) 2408 zz=zz+wtatcon(ii,iatom)**2/mass(iatom) 2409 end do 2410 end do 2411 end do 2412 vec(:)=force_centroid(:)-weightsum(:)*af/zz 2413 do ii=1,3 2414 mat(:,ii)=-weightsum(:)*weightsum(ii)/zz 2415 mat(ii,ii)=mat(ii,ii)+masstot 2416 end do 2417 call matr3inv(mat,matinv) 2418 2419 !Calculation of a Lagrange multipliers: 2420 ! lambda_cst: to apply the constraint 2421 ! lambda_com: to maintain the position of the center of mass 2422 lambda_com(:)=matmul(matinv,vec) 2423 lambda_cst=(af-dot_product(weightsum,lambda_com))*dble(trotter)/zz 2424 2425 !Modification of forces 2426 one_over_trotter=one/dble(trotter) 2427 do iimage=1,trotter 2428 do iatom=1,natom 2429 do ii=1,3 2430 forces(ii,iatom,iimage)=forces(ii,iatom,iimage) & 2431 & -lambda_cst*wtatcon(ii,iatom)*one_over_trotter & 2432 & -lambda_com(ii)*mass(iatom) 2433 end do 2434 end do 2435 end do 2436 2437 !Computation of relevant outputs 2438 constraint_output(:)=zero 2439 !1-Reaction coordinate 2440 do iatom=1,natom 2441 do ii=1,3 2442 xcart_centroid=zero 2443 do iimage=1,trotter 2444 xcart_centroid=xcart_centroid+xcart(ii,iatom,iimage) 2445 end do 2446 constraint_output(1)=constraint_output(1)+xcart_centroid*wtatcon(ii,iatom) 2447 end do 2448 end do 2449 constraint_output(1)=constraint_output(1)/dble(trotter) 2450 !2-Force on reaction coordinate 2451 constraint_output(2)=-lambda_cst 2452 2453 end select 2454 2455 end subroutine pimd_apply_constraint
m_pimd/pimd_coord_transform [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_coord_transform
FUNCTION
Apply a coordinate transformation on a given vector field (defined for each atom in in cell) - for PIMD Possible choices for the transformation: 0: no transformation 1: normal mode tranformation 2: staging transformation
INPUTS
ioption=option given the direction of the transformation +1: from primitive coordinates to transformed coordinates -1: from transformed coordinates to primitive coordinates natom=number of atoms transform=coordinate transformation: 0: no tranformation 1: normal mode transformation 2: staging transformation trotter=Trotter number
OUTPUT
SIDE EFFECTS
array(3,natom,trotter)=array to be transformed
SOURCE
2064 subroutine pimd_coord_transform(array,ioption,natom,transform,trotter) 2065 2066 !Arguments ------------------------------------ 2067 !scalars 2068 integer,intent(in) :: ioption,natom,transform,trotter 2069 !arrays 2070 real(dp),intent(inout) :: array(3,natom,trotter) 2071 !Local variables------------------------------- 2072 !scalars 2073 integer :: iatom,ii,iimage,iimagem,iimagep,iimage2,ll 2074 !arrays 2075 real(dp),allocatable :: array_temp(:,:,:) 2076 real(dp),allocatable :: nrm(:,:) 2077 2078 !************************************************************************ 2079 2080 !=== No transformation =================================================== 2081 if (transform==0) then 2082 2083 return 2084 2085 !=== Normal mode transformation ========================================== 2086 else if (transform==1) then 2087 2088 ! ---------------- From primitive to transformed coordinates ------------ 2089 if (ioption==+1) then 2090 ABI_MALLOC(array_temp,(3,natom,trotter)) 2091 2092 array_temp(:,:,1)=zero 2093 do iimage=1,trotter 2094 array_temp(:,:,1)=array_temp(:,:,1)+array(:,:,iimage) 2095 end do 2096 2097 do ll=2,trotter/2 2098 array_temp(:,:,2*ll-2)=zero 2099 array_temp(:,:,2*ll-1)=zero 2100 do iimage=1,trotter 2101 array_temp(:,:,2*ll-2)=array_temp(:,:,2*ll-2)+array(:,:,iimage)* & 2102 & cos(two_pi*dble((iimage-1)*(ll-1))/dble(trotter)) 2103 array_temp(:,:,2*ll-1)=array_temp(:,:,2*ll-1)-array(:,:,iimage)* & 2104 & sin(two_pi*dble((iimage-1)*(ll-1))/dble(trotter)) 2105 end do 2106 end do 2107 2108 array_temp(:,:,trotter)=zero 2109 do iimage=1,trotter 2110 array_temp(:,:,trotter)=array_temp(:,:,trotter)+& 2111 & array(:,:,iimage)*dble((-1)**(iimage-1)) 2112 end do 2113 2114 array=array_temp/dble(trotter) 2115 2116 ABI_FREE(array_temp) 2117 2118 ! ---------------- From transformed to primitive coordinates ------------ 2119 else if (ioption==-1) then 2120 2121 ABI_MALLOC(array_temp,(3,natom,trotter)) !real part 2122 ABI_MALLOC(nrm,(trotter,trotter)) !real part 2123 2124 do iimage=1,trotter 2125 nrm(iimage,1)=one 2126 nrm(iimage,trotter)=dble((-1)**(iimage-1)) 2127 end do 2128 2129 do ll=2,trotter/2 2130 do iimage=1,trotter 2131 nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* & 2132 & (iimage-1))/dble(trotter)) 2133 nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* & 2134 & (iimage-1))/dble(trotter)) 2135 end do 2136 end do 2137 2138 do iimage=1,trotter 2139 array_temp(:,:,iimage)=zero 2140 do ll=1,trotter 2141 array_temp(:,:,iimage)=array_temp(:,:,iimage)+nrm(iimage,ll)*array(:,:,ll) 2142 end do 2143 end do 2144 array=array_temp 2145 2146 ABI_FREE(array_temp) 2147 ABI_FREE(nrm) 2148 2149 end if ! ioption 2150 2151 !=== Staging transformation ============================================== 2152 else if (transform==2) then 2153 2154 ! ---------------- From primitive to transformed coordinates ------------ 2155 if (ioption==+1) then 2156 2157 ABI_MALLOC(array_temp,(3,natom,trotter)) 2158 array_temp=zero 2159 do iimage=1,trotter 2160 iimagep=iimage+1;if(iimage==trotter) iimagep=1 2161 iimagem=iimage-1;if(iimage==1) iimagem=trotter 2162 do iatom=1,natom 2163 do ii=1,3 2164 array_temp(ii,iatom,iimage)=(dble(iimagem)*array(ii,iatom,iimagep)+array(ii,iatom,1))/dble(iimage) 2165 end do 2166 end do 2167 end do 2168 if (trotter>1) then 2169 do iimage=2,trotter 2170 do iatom=1,natom 2171 do ii=1,3 2172 array(ii,iatom,iimage)=array(ii,iatom,iimage)-array_temp(ii,iatom,iimage) 2173 end do 2174 end do 2175 end do 2176 end if 2177 ABI_FREE(array_temp) 2178 2179 ! ---------------- From transformed to primitive coordinates ------------ 2180 else if (ioption==-1) then 2181 2182 ABI_MALLOC(array_temp,(3,natom,trotter)) 2183 array_temp=zero 2184 do iimage=1,trotter 2185 do iatom=1,natom 2186 do ii=1,3 2187 array_temp(ii,iatom,iimage)=array(ii,iatom,1) 2188 end do 2189 end do 2190 end do 2191 if (trotter>1) then 2192 do iimage=2,trotter 2193 do iatom=1,natom 2194 do ii=1,3 2195 do iimage2=iimage,trotter 2196 array_temp(ii,iatom,iimage)=array_temp(ii,iatom,iimage) & 2197 & +array(ii,iatom,iimage2)*(dble(iimage-1))/(dble(iimage2-1)) 2198 end do 2199 end do 2200 end do 2201 end do 2202 end if 2203 array(:,:,:)=array_temp(:,:,:) 2204 ABI_FREE(array_temp) 2205 2206 end if ! ioption 2207 2208 end if ! transform 2209 2210 end subroutine pimd_coord_transform
m_pimd/pimd_destroy [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_destroy
FUNCTION
Destroy the content of a datastructure of type pimd_type. Close open file(s) related to this datastructure.
INPUTS
OUTPUT
SIDE EFFECTS
pimd_param=datastructure of type pimd_type. several parameters for PIMD.
SOURCE
275 subroutine pimd_destroy(pimd_param) 276 277 implicit none 278 279 !Arguments ------------------------------------ 280 !scalars 281 type(pimd_type),intent(inout) :: pimd_param 282 !arrays 283 !Local variables------------------------------- 284 !scalars 285 integer :: ierr 286 character(len=100) :: msg 287 !arrays 288 289 !************************************************************************ 290 291 ABI_SFREE(pimd_param%zeta_prev) 292 ABI_SFREE(pimd_param%zeta) 293 ABI_SFREE(pimd_param%zeta_next) 294 ABI_SFREE(pimd_param%dzeta) 295 296 if (pimd_param%qtb_file_unit>0) then 297 if (is_open(pimd_param%qtb_file_unit)) then 298 ierr=close_unit(pimd_param%qtb_file_unit,msg) 299 end if 300 end if 301 302 if (pimd_param%traj_unit>0) then 303 if (is_open(pimd_param%traj_unit)) then 304 ierr=close_unit(pimd_param%traj_unit,msg) 305 end if 306 end if 307 308 call pimd_nullify(pimd_param) 309 310 end subroutine pimd_destroy
m_pimd/pimd_diff_stress [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_diff_stress
FUNCTION
Compute the difference between the stress tensor and the stress target
INPUTS
stress_pimd(3,3,3)=stress tensor for PIMD Last dimension (3) corresponds to 3 different pressure estimators stress_target(6)=stress target
OUTPUT
pimd_diff_stress(3,3,3)=difference between stresses and stress target First dimension (3) corresponds to 3 different pressure estimators
SIDE EFFECTS
SOURCE
1696 function pimd_diff_stress(stress_pimd,stress_target) 1697 1698 !Arguments ------------------------------------ 1699 !scalars 1700 !arrays 1701 real(dp),intent(in) :: stress_pimd(3,3,3),stress_target(6) 1702 real(dp) :: pimd_diff_stress(3,3) 1703 !Local variables------------------------------- 1704 !scalars 1705 integer :: ii,jj 1706 !arrays 1707 real(dp) :: stress_pimd2(3,3) 1708 1709 !************************************************************************ 1710 1711 !Choice: the primitive estimator for pressure is chosen 1712 do ii=1,3 1713 do jj=1,3 1714 stress_pimd2(ii,jj)=stress_pimd(1,ii,jj) 1715 end do 1716 end do 1717 1718 !+stress_target instead of - because it is translated from stress to pressure tensor 1719 pimd_diff_stress(1,1)=stress_pimd2(1,1)+stress_target(1) 1720 pimd_diff_stress(2,2)=stress_pimd2(2,2)+stress_target(2) 1721 pimd_diff_stress(3,3)=stress_pimd2(3,3)+stress_target(3) 1722 pimd_diff_stress(2,3)=stress_pimd2(2,3)+stress_target(4) 1723 pimd_diff_stress(3,2)=stress_pimd2(3,2)+stress_target(4) 1724 pimd_diff_stress(1,3)=stress_pimd2(1,3)+stress_target(5) 1725 pimd_diff_stress(3,1)=stress_pimd2(3,1)+stress_target(5) 1726 pimd_diff_stress(1,2)=stress_pimd2(1,2)+stress_target(6) 1727 pimd_diff_stress(2,1)=stress_pimd2(2,1)+stress_target(6) 1728 1729 end function pimd_diff_stress
m_pimd/pimd_energies [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_energies
FUNCTION
In the case od PIMD, compute the several contribution to total energy
INPUTS
etotal_img(trotter)= energy (from DFT) for each cell forces(3,natom,trotter)=forces (from DFT) on atoms in each cell natom=number of atoms spring(natom)=spring constants in the primitive scheme trotter=Trotter number xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
OUTPUT
eharm =harmonic energy eharm_virial=harmonic energy from virial estimator epot =potential energy
SIDE EFFECTS
SOURCE
1260 subroutine pimd_energies(eharm,eharm_virial,epot,etotal_img,forces,natom,spring,trotter,xcart) 1261 1262 !Arguments ------------------------------------ 1263 !scalars 1264 integer,intent(in) :: natom,trotter 1265 real(dp),intent(out) :: eharm,eharm_virial,epot 1266 !arrays 1267 real(dp),intent(in) :: etotal_img(trotter),forces(3,natom,trotter) 1268 real(dp),intent(in) :: xcart(3,natom,trotter) 1269 real(dp),intent(in) :: spring(natom) 1270 !Local variables------------------------------- 1271 !scalars 1272 integer :: iatom,ii,iimage,iimagep 1273 !arrays 1274 real(dp),allocatable :: centroid(:,:) 1275 1276 !************************************************************************ 1277 1278 !Compute the centroid 1279 ABI_MALLOC(centroid,(3,natom)) 1280 centroid=zero 1281 do iimage=1,trotter 1282 do iatom=1,natom 1283 do ii=1,3 1284 centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage) 1285 end do 1286 end do 1287 end do 1288 centroid=centroid/dble(trotter) 1289 1290 !Potential energy 1291 epot=sum(etotal_img(1:trotter))/dble(trotter) 1292 1293 !Harmonic energy 1294 eharm=zero 1295 do iimage=1,trotter 1296 iimagep=iimage+1;if(iimage==trotter)iimagep=1 1297 do iatom=1,natom 1298 do ii=1,3 1299 eharm=eharm+half*spring(iatom)*((xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))**2) 1300 end do 1301 end do 1302 end do 1303 1304 !Harmonic energy from virial estimator 1305 eharm_virial=zero 1306 do iimage=1,trotter 1307 do iatom=1,natom 1308 do ii=1,3 1309 eharm_virial=eharm_virial-(xcart(ii,iatom,iimage)-centroid(ii,iatom)) & 1310 & *forces(ii,iatom,iimage) 1311 end do 1312 end do 1313 end do 1314 eharm_virial=eharm_virial/dble(two*trotter) 1315 1316 ABI_FREE(centroid) 1317 1318 end subroutine pimd_energies
m_pimd/pimd_force_transform [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_force_transform
FUNCTION
Apply a coordinate transformation on forces (defined for each atom in in cell) - for PIMD Possible choices for the transformation: 0: no transformation 1: normal mode tranformation 2: staging transformation
INPUTS
ioption=option given the direction of the transformation +1: from primitive coordinates to transformed coordinates -1: from transformed coordinates to primitive coordinates natom=number of atoms transform=coordinate transformation: 0: no tranformation 1: normal mode transformation 2: staging transformation trotter=Trotter number
OUTPUT
SIDE EFFECTS
forces(3,natom,trotter)=array containing forces
NOTES
Back transformation (ioption=-1) not implemented !
SOURCE
2247 subroutine pimd_force_transform(forces,ioption,natom,transform,trotter) 2248 2249 !Arguments ------------------------------------ 2250 !scalars 2251 integer,intent(in) :: ioption,natom,transform,trotter 2252 !arrays 2253 real(dp),intent(inout) :: forces(3,natom,trotter) 2254 !Local variables------------------------------- 2255 !scalars 2256 integer :: iatom,ii,iimage,ll 2257 character(len=500) :: msg 2258 !arrays 2259 real(dp),allocatable :: forces_temp(:,:,:),nrm(:,:) 2260 2261 !************************************************************************ 2262 2263 if (ioption==-1) then 2264 msg='Back transformation not implemented !' 2265 ABI_BUG(msg) 2266 end if 2267 2268 !=== No transformation =================================================== 2269 if (transform==0) then 2270 2271 return 2272 2273 !=== Normal mode transformation ========================================== 2274 else if (transform==1) then 2275 2276 !normal mode forces 2277 ABI_MALLOC(forces_temp,(3,natom,trotter)) 2278 ABI_MALLOC(nrm,(trotter,trotter)) 2279 2280 do iimage=1,trotter 2281 nrm(iimage,1)=one 2282 nrm(iimage,trotter)=dble((-1)**(iimage-1)) 2283 end do 2284 2285 do ll=2,trotter/2 2286 do iimage=1,trotter 2287 nrm(iimage,2*ll-2)= two*cos(two_pi*dble((ll-1)* & 2288 & (iimage-1))/dble(trotter)) 2289 nrm(iimage,2*ll-1)=-two*sin(two_pi*dble((ll-1)* & 2290 & (iimage-1))/dble(trotter)) 2291 end do 2292 end do 2293 2294 do ll=1,trotter 2295 forces_temp(:,:,ll)=zero 2296 do iimage=1,trotter 2297 forces_temp(:,:,ll)=forces_temp(:,:,ll)+nrm(iimage,ll)*forces(:,:,iimage) 2298 end do 2299 end do 2300 2301 forces=forces_temp 2302 2303 ABI_FREE(forces_temp) 2304 ABI_FREE(nrm) 2305 2306 !=== Staging transformation ============================================== 2307 else if (transform==2) then 2308 2309 !staging forces 2310 ABI_MALLOC(forces_temp,(3,natom,trotter)) 2311 forces_temp=zero 2312 do iimage=1,trotter 2313 do iatom=1,natom 2314 do ii=1,3 2315 forces_temp(ii,iatom,1)=forces_temp(ii,iatom,1)+forces(ii,iatom,iimage) 2316 end do 2317 end do 2318 end do 2319 if (trotter>1) then 2320 do iimage=2,trotter 2321 do iatom=1,natom 2322 do ii=1,3 2323 forces_temp(ii,iatom,iimage)=forces(ii,iatom,iimage) & 2324 & +forces_temp(ii,iatom,iimage-1)*(dble(iimage-2)/dble(iimage-1)) 2325 end do 2326 end do 2327 end do 2328 end if 2329 forces=forces_temp 2330 ABI_FREE(forces_temp) 2331 2332 end if ! transform 2333 2334 end subroutine pimd_force_transform
m_pimd/pimd_forces [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_forces
FUNCTION
Modify forces in order to take into account PIMD contribution
INPUTS
natom=number of atoms spring(natom,spring_dim)=spring constants (spring_dim=1 or trotter) transform=coordinate transformation: 0: no tranformation 1: normal mode transformation 2: staging transformation trotter=Trotter number xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell
OUTPUT
SIDE EFFECTS
forces(3,natom,trotter)= at input: forces from electronic calculation at output: forces from electronic calculation + quantum spring contribution
SOURCE
1349 subroutine pimd_forces(forces,natom,spring,transform,trotter,xcart) 1350 1351 !Arguments ------------------------------------ 1352 !scalars 1353 integer,intent(in) :: natom,transform,trotter 1354 !arrays 1355 real(dp),intent(in) :: xcart(3,natom,trotter) 1356 real(dp),intent(in) :: spring(:,:) 1357 real(dp),intent(inout) :: forces(3,natom,trotter) 1358 !Local variables------------------------------- 1359 !scalars 1360 integer :: iatom,ii,iimage,iimagem,iimagep,ispring,natom_spring,nspring 1361 character(len=500) :: msg 1362 !arrays 1363 1364 !************************************************************************ 1365 1366 natom_spring=size(spring,1);nspring=size(spring,2) 1367 if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then 1368 msg='Wrong dimensions for array spring !' 1369 ABI_BUG(msg) 1370 end if 1371 1372 if (transform==0) then 1373 do iimage=1,trotter 1374 ispring=min(nspring,iimage) 1375 iimagep=iimage+1; iimagem=iimage-1 1376 if(iimage==trotter) iimagep=1 1377 if(iimage==1) iimagem=trotter 1378 do iatom=1,natom 1379 do ii=1,3 1380 forces(ii,iatom,iimage)= & 1381 & forces(ii,iatom,iimage)/dble(trotter) & 1382 & - spring(iatom,ispring)*(two*xcart(ii,iatom,iimage)-xcart(ii,iatom,iimagem) & 1383 & -xcart(ii,iatom,iimagep)) 1384 end do 1385 end do 1386 end do 1387 1388 else 1389 do iimage=1,trotter 1390 ispring=min(nspring,iimage) 1391 do iatom=1,natom 1392 do ii=1,3 1393 forces(ii,iatom,iimage)= & 1394 & forces(ii,iatom,iimage)/dble(trotter) & 1395 & - spring(iatom,ispring)*xcart(ii,iatom,iimage) 1396 end do 1397 end do 1398 end do 1399 1400 end if 1401 1402 end subroutine pimd_forces
m_pimd/pimd_init [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_init
FUNCTION
Initialize a datastructure of type pimd_type. Open file(s) related to this datastructure.
INPUTS
dtset <type(dataset_type)>=all input variables in current dataset is_master= TRUE if I am the master process (proc 0)
OUTPUT
SIDE EFFECTS
pimd_param=datastructure of type pimd_type. several parameters for Path-Integral MD.
SOURCE
135 subroutine pimd_init(dtset,pimd_param,is_master) 136 137 implicit none 138 139 !Arguments ------------------------------------ 140 !scalars 141 logical,intent(in) :: is_master 142 type(dataset_type),target,intent(in) :: dtset 143 type(pimd_type),intent(inout) :: pimd_param 144 !Local variables------------------------------- 145 !scalars 146 integer :: ierr 147 character(len=200) :: msg 148 149 !************************************************************************ 150 151 call pimd_nullify(pimd_param) 152 153 if((dtset%imgmov==9).or.(dtset%imgmov==10).or.(dtset%imgmov==13))then 154 pimd_param%adpimd = dtset%adpimd 155 pimd_param%constraint = dtset%pimd_constraint 156 pimd_param%irandom = dtset%irandom 157 pimd_param%nnos = dtset%nnos 158 pimd_param%ntypat = dtset%ntypat 159 pimd_param%optcell = dtset%optcell 160 pimd_param%pitransform = dtset%pitransform 161 pimd_param%adpimd_gamma= dtset%adpimd_gamma 162 pimd_param%vis = dtset%vis 163 pimd_param%bmass = dtset%bmass 164 pimd_param%dtion = dtset%dtion 165 pimd_param%friction = dtset%friction 166 pimd_param%mdtemp =>dtset%mdtemp 167 pimd_param%pimass =>dtset%pimass 168 pimd_param%strtarget =>dtset%strtarget 169 pimd_param%amu =>dtset%amu_orig(:,1) 170 pimd_param%qmass =>dtset%qmass 171 pimd_param%typat =>dtset%typat 172 pimd_param%wtatcon =>dtset%wtatcon 173 if(dtset%imgmov==10)then 174 pimd_param%use_qtb=1 175 if(is_master)then 176 call pimd_init_qtb(dtset,pimd_param%qtb_file_unit) 177 end if 178 end if 179 if(dtset%imgmov==13)then 180 ABI_MALLOC(pimd_param%zeta_prev,(3,dtset%natom,dtset%nimage,dtset%nnos)) 181 ABI_MALLOC(pimd_param%zeta ,(3,dtset%natom,dtset%nimage,dtset%nnos)) 182 ABI_MALLOC(pimd_param%zeta_next,(3,dtset%natom,dtset%nimage,dtset%nnos)) 183 ABI_MALLOC(pimd_param%dzeta ,(3,dtset%natom,dtset%nimage,dtset%nnos)) 184 pimd_param%zeta_prev=zero 185 pimd_param%zeta =zero 186 pimd_param%zeta_next=zero 187 pimd_param%dzeta =zero 188 end if 189 if(dtset%useria==37)then 190 ierr=open_file('pimd_traj.dat',msg,newunit=pimd_param%traj_unit,form='unformatted') 191 if (ierr/=0) then 192 ABI_ERROR(msg) 193 end if 194 end if 195 end if 196 197 end subroutine pimd_init
m_pimd/pimd_init_qtb [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_init_qtb
FUNCTION
Only relevant for PIMD + Quantum Thermal Bath (QTB); Initialize reading of PIQTB random force file. This routine should be called only by master proc.
INPUTS
dtset <type(dataset_type)>=all input variables in current dataset
OUTPUT
qtb_file_unit=if a PIQTB_force file exists, return its file unit.
SOURCE
332 subroutine pimd_init_qtb(dtset,qtb_file_unit) 333 334 implicit none 335 336 !Arguments ------------------------------------ 337 !scalars 338 integer,intent(out) :: qtb_file_unit 339 type(dataset_type),target,intent(in) :: dtset 340 !Local variables------------------------------- 341 !scalars 342 integer :: ierr,ndof_qtb,ntimimage_qtb,nimage_qtb 343 real(dp) :: dtion_qtb,mdtemp_qtb 344 character(len=200) :: msg 345 346 !************************************************************************ 347 348 !Try to open PIQTB random force file 349 ierr=open_file('piqtb_force',msg,newunit=qtb_file_unit,& 350 & form='unformatted',status='old') 351 352 !Read first line of the file 353 read(qtb_file_unit) dtion_qtb,ntimimage_qtb,mdtemp_qtb,nimage_qtb,ndof_qtb 354 355 !Check consistency of the read parameters with ABINIT input file 356 if (abs(dtion_qtb-dtset%dtion)>tol6) then 357 msg='dtion read from piqtb_force file different from dtion in input file!' 358 ABI_ERROR(msg) 359 end if 360 if (abs(mdtemp_qtb-dtset%mdtemp(2))>tol6) then 361 msg='mdtemp read from piqtb_force file different from mdtemp(2) in input file!' 362 ABI_ERROR(msg) 363 end if 364 if (ntimimage_qtb<dtset%ntimimage) then 365 msg='ntimimage read from piqtb_force file smaller than ntimimage in input file!' 366 ABI_ERROR(msg) 367 end if 368 if (nimage_qtb/=dtset%nimage) then 369 msg='nimage read from piqtb_force file different from nimage in input file!' 370 ABI_ERROR(msg) 371 end if 372 if (ndof_qtb/=3*dtset%natom*dtset%nimage) then 373 msg='Nb of degrees of freedom read from piqtb_force not consistent with input file!' 374 ABI_ERROR(msg) 375 end if 376 377 end subroutine pimd_init_qtb
m_pimd/pimd_initvel [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_initvel
FUNCTION
Initialize velocities for PIMD with a gaussian distribution fixing the center of mass and eventually applying a constraint on atomic positions
INPUTS
constraint=type of constraint to be applied mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms temperature=temperature used to define velocities trotter=Trotter number wtatcon(3,natom)=weights for atomic constraints
OUTPUT
vel(3,natom,trotter)=velocities of atoms in in each cell
SIDE EFFECTS
iseed=seed for random number generator
SOURCE
838 subroutine pimd_initvel(iseed,mass,natom,temperature,trotter,vel,constraint,wtatcon) 839 840 !Arguments ------------------------------------ 841 !scalars 842 integer,intent(in) :: constraint,natom,trotter 843 integer,intent(inout) :: iseed 844 real(dp),intent(in) :: temperature 845 !arrays 846 real(dp),intent(in) :: mass(:,:),wtatcon(3,natom) 847 real(dp),intent(out) :: vel(3,natom,trotter) 848 !Local variables------------------------------- 849 !scalars 850 integer :: iatom,ii,iimage,imass,natom_mass,nmass 851 real(dp) :: mtot,rescale_vel 852 character(len=500) :: msg 853 !arrays 854 real(dp) :: mvini(3) 855 856 !************************************************************************ 857 858 natom_mass=size(mass,1);nmass=size(mass,2) 859 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 860 msg='Wrong dimensions for array mass !' 861 ABI_BUG(msg) 862 end if 863 864 !Compute total mass (of non-constrained atoms) 865 if (constraint==0) then 866 mtot=sum(mass(1:natom,1:nmass)) 867 else 868 mtot=zero 869 do iatom=1,natom 870 if (all(abs(wtatcon(:,iatom))<tol8)) mtot=mtot+sum(mass(iatom,1:nmass)) 871 end do 872 end if 873 if (nmass==1) mtot=mtot*dble(trotter) 874 875 !Initialize randomly the velocities 876 do iimage=1,trotter 877 imass=min(nmass,iimage) 878 do iatom=1,natom 879 do ii=1,3 880 vel(ii,iatom,iimage)=sqrt(kb_HaK*temperature/mass(iatom,imass))*cos(two_pi*uniformrandom(iseed)) 881 vel(ii,iatom,iimage)=vel(ii,iatom,iimage)*sqrt(-two*log(uniformrandom(iseed))) 882 end do 883 end do 884 end do 885 886 !Cancel velocities of constrained atoms 887 if (constraint/=0) then 888 do iimage=1,trotter 889 do iatom=1,natom 890 if (any(abs(wtatcon(:,iatom))>=tol8)) vel(:,iatom,iimage)=zero 891 end do 892 end do 893 end if 894 895 !Make sure that the (sum of m_i v_i) at step zero is zero 896 mvini=zero 897 do iimage=1,trotter 898 imass=min(nmass,iimage) 899 do iatom=1,natom 900 do ii=1,3 901 mvini(ii)=mvini(ii)+mass(iatom,imass)*vel(ii,iatom,iimage) 902 end do 903 end do 904 end do 905 if (constraint==0) then 906 do iimage=1,trotter 907 do iatom=1,natom 908 do ii=1,3 909 vel(ii,iatom,iimage)=vel(ii,iatom,iimage)-(mvini(ii)/mtot) 910 end do 911 end do 912 end do 913 else 914 do iimage=1,trotter 915 do iatom=1,natom 916 if (all(abs(wtatcon(:,iatom))<tol8)) vel(:,iatom,iimage)=vel(:,iatom,iimage)-(mvini(:)/mtot) 917 end do 918 end do 919 end if 920 921 !Now rescale the velocities to give the exact temperature 922 rescale_vel=sqrt(temperature/pimd_temperature(mass,vel)) 923 vel(:,:,:)=vel(:,:,:)*rescale_vel 924 925 end subroutine pimd_initvel
m_pimd/pimd_is_restart [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_is_restart
FUNCTION
Determine whether this is a PIMD restart or not: test on value of velocities and corresponding temperature
INPUTS
mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) vel(3,natom,nimage)=velocities for each image of the cell
OUTPUT
pimd_is_restart=1 if temperature is not zero
SIDE EFFECTS
SOURCE
445 function pimd_is_restart(mass,vel,vel_cell) 446 447 !Arguments ------------------------------------ 448 !scalars 449 integer :: pimd_is_restart 450 !arrays 451 real(dp),intent(in) :: mass(:,:),vel(:,:,:) 452 real(dp),intent(in),optional :: vel_cell(:,:) 453 !Local variables------------------------------- 454 !scalars 455 real(dp),parameter :: zero_temp=tol7 456 !arrays 457 458 !************************************************************************ 459 460 pimd_is_restart=0 461 if (pimd_temperature(mass,vel)>zero_temp) pimd_is_restart=1 462 if (present(vel_cell)) then 463 if (maxval(vel_cell)>zero_temp) pimd_is_restart=pimd_is_restart+10 464 end if 465 466 end function pimd_is_restart
m_pimd/pimd_langevin_forces [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_langevin_forces
FUNCTION
Compute Langevin contribution to PIMD forces
INPUTS
alea(3,natom,trotter)=set of random numbers forces(3,natom,trotter)=forces without Langevin contribution friction=friction factor langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter) mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms trotter=Trotter number vel(3,natom,trotter)=velocities of atoms in each cell
OUTPUT
forces_langevin(3,natom,trotter)=forces including Langevin contribution
SIDE EFFECTS
SOURCE
1431 subroutine pimd_langevin_forces(alea,forces,forces_langevin,friction,& 1432 & langev,mass,natom,trotter,vel) 1433 1434 !Arguments ------------------------------------ 1435 !scalars 1436 integer,intent(in) :: natom,trotter 1437 real(dp),intent(in) :: friction 1438 !arrays 1439 real(dp),intent(in) :: alea(3,natom,trotter),forces(3,natom,trotter) 1440 real(dp),intent(in) :: langev(:,:),mass(:,:),vel(3,natom,trotter) 1441 real(dp),intent(out) :: forces_langevin(3,natom,trotter) 1442 !Local variables------------------------------- 1443 !scalars 1444 integer :: iatom,ii,iimage,imass,natom_mass,nmass 1445 character(len=500) :: msg 1446 !arrays 1447 1448 !************************************************************************ 1449 1450 natom_mass=size(mass,1);nmass=size(mass,2) 1451 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1452 msg='Wrong dimensions for array mass !' 1453 ABI_BUG(msg) 1454 end if 1455 1456 do iimage=1,trotter 1457 imass=min(nmass,iimage) 1458 do iatom=1,natom 1459 do ii=1,3 1460 forces_langevin(ii,iatom,iimage)=forces(ii,iatom,iimage) & 1461 & + langev(iatom,imass)*alea(ii,iatom,iimage) & 1462 & - friction*mass(iatom,imass)*vel(ii,iatom,iimage) 1463 end do 1464 end do 1465 end do 1466 1467 end subroutine pimd_langevin_forces
m_pimd/pimd_langevin_random [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_langevin_random
FUNCTION
Generate a set of random numbers to be used for PIMD Langevin algorithm
INPUTS
irandom=option for random number generator: 1:uniform random routine provided within Abinit package 2:Fortran 90 random number generator 3:ZBQLU01 non deterministic random number generator langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter) mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms trotter=Trotter number zeroforce=flag; if 1 keep sum of forces equal to zero
OUTPUT
alea(3,natom,trotter)=set of random numbers
SIDE EFFECTS
iseed=seed for random number generator (used only if irandom=1)
SOURCE
956 subroutine pimd_langevin_random(alea,irandom,iseed,langev,mass,natom,trotter,zeroforce) 957 958 !Arguments ------------------------------------ 959 !scalars 960 integer,intent(in) :: irandom,natom,trotter,zeroforce 961 integer,intent(inout) :: iseed 962 !arrays 963 real(dp),intent(in) :: langev(:,:),mass(:,:) 964 real(dp),intent(out) :: alea(3,natom,trotter) 965 !Local variables------------------------------- 966 !scalars 967 integer :: iatom,ii,iimage,imass,nmass,natom_mass 968 real(dp) :: mtot,r1,r2 969 character(len=500) :: msg 970 !arrays 971 real(dp) :: total(3) 972 973 !************************************************************************ 974 975 natom_mass=size(mass,1);nmass=size(mass,2) 976 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 977 msg='Wrong dimensions for array mass !' 978 ABI_BUG(msg) 979 end if 980 981 mtot=sum(mass(1:natom,1:nmass)) 982 if (nmass==1) mtot=mtot*trotter 983 984 !Draw random numbers 985 do iimage=1,trotter 986 do iatom=1,natom 987 do ii=1,3 988 select case(irandom) 989 case(1) 990 r1=uniformrandom(iseed) 991 r2=uniformrandom(iseed) 992 case(2) 993 call random_number(r1) 994 call random_number(r2) 995 case(3) 996 r1=ZBQLU01(zero) 997 r2=ZBQLU01(zero) 998 end select 999 alea(ii,iatom,iimage)= cos(two_pi*r1)*sqrt(-log(r2)*two) 1000 end do 1001 end do 1002 end do 1003 1004 !Make sure that the sum of random forces is zero 1005 if(zeroforce==1)then 1006 total=zero 1007 mtot=zero 1008 do iimage=1,trotter 1009 imass=min(nmass,iimage) 1010 do iatom=1,natom 1011 mtot=mtot+mass(iatom,imass) 1012 end do 1013 end do 1014 do iimage=1,trotter 1015 imass=min(nmass,iimage) 1016 do iatom=1,natom 1017 do ii=1,3 1018 total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage) 1019 end do 1020 end do 1021 end do 1022 do iimage=1,trotter 1023 imass=min(nmass,iimage) 1024 do iatom=1,natom 1025 do ii=1,3 1026 alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- & 1027 & (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot) 1028 end do 1029 end do 1030 end do 1031 ! now random forces have been rescaled so that their sum is zero 1032 end if 1033 1034 end subroutine pimd_langevin_random
m_pimd/pimd_langevin_random_bar [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_langevin_random_bar
FUNCTION
Generate a set of random numbers to be used for the barostat of PIMD Langevin algorithm
INPUTS
irandom=option for random number generator: 1:uniform random routine provided within Abinit package 2:Fortran 90 random number generator 3:ZBQLU01 non deterministic random number generator
OUTPUT
alea_bar(3,3)=set of random numbers
SIDE EFFECTS
iseed=seed for random number generator (used only if irandom=1)
SOURCE
1148 subroutine pimd_langevin_random_bar(alea_bar,irandom,iseed) 1149 1150 !Arguments ------------------------------------ 1151 !scalars 1152 integer,intent(in) :: irandom 1153 integer,intent(inout) :: iseed 1154 !arrays 1155 real(dp),intent(out) :: alea_bar(3,3) 1156 !Local variables------------------------------- 1157 !scalars 1158 integer :: ii,jj 1159 real(dp) :: r1,r2 1160 !arrays 1161 1162 !************************************************************************ 1163 1164 !Draw random numbers 1165 do ii=1,3 1166 do jj=1,3 1167 select case(irandom) 1168 case(1) 1169 r1=uniformrandom(iseed) 1170 r2=uniformrandom(iseed) 1171 case(2) 1172 call random_number(r1) 1173 call random_number(r2) 1174 case(3) 1175 r1=ZBQLU01(zero) 1176 r2=ZBQLU01(zero) 1177 end select 1178 alea_bar(ii,jj)= cos(two_pi*r1)*sqrt(-log(r2)*two) 1179 end do 1180 end do 1181 1182 !Symmetrize 1183 alea_bar(1,2)=half*(alea_bar(1,2)+alea_bar(2,1)) 1184 alea_bar(1,3)=half*(alea_bar(1,3)+alea_bar(3,1)) 1185 alea_bar(2,3)=half*(alea_bar(2,3)+alea_bar(3,2)) 1186 alea_bar(2,1)=alea_bar(1,2) 1187 alea_bar(3,1)=alea_bar(1,3) 1188 alea_bar(3,2)=alea_bar(2,3) 1189 1190 end subroutine pimd_langevin_random_bar
m_pimd/pimd_langevin_random_init [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_langevin_random_init
FUNCTION
Initialize random number generator to be used for PIMD Langevin algorithm
INPUTS
irandom=option for random number generator: 1:uniform random routine provided within Abinit package 2:Fortran 90 random number generator 3:ZBQLU01 non deterministic random number generator
OUTPUT
SIDE EFFECTS
iseed=seed for random number generator (used only if irandom=1)
SOURCE
1215 subroutine pimd_langevin_random_init(irandom,iseed) 1216 1217 !Arguments ------------------------------------ 1218 !scalars 1219 integer,intent(in) :: irandom 1220 integer,intent(inout) :: iseed 1221 1222 !************************************************************************ 1223 1224 if (irandom==3) then 1225 call ZBQLINI(0) 1226 end if 1227 1228 !Fake statement 1229 return;if (.false.) iseed=zero 1230 1231 end subroutine pimd_langevin_random_init
m_pimd/pimd_langevin_random_qtb [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_langevin_random_qtb
FUNCTION
Read a set of random forces (atm units) to be used for PIMD QTB algorithm
INPUTS
langev(natom,mass_dim)=Langevin factors (mass_dim=1 or trotter) mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms qtb_file_unit= ramdom forces file unit trotter=Trotter number zeroforce=flag; if 1 keep sum of forces equal to zero
OUTPUT
alea(3,natom,trotter)=set of random forces
SOURCE
1059 subroutine pimd_langevin_random_qtb(alea,langev,mass,natom,qtb_file_unit,trotter,zeroforce) 1060 1061 !Arguments ------------------------------------ 1062 !scalars 1063 integer,intent(in) :: natom,qtb_file_unit,trotter,zeroforce 1064 !arrays 1065 real(dp),intent(in) :: langev(:,:),mass(:,:) 1066 real(dp),intent(out) :: alea(3,natom,trotter) 1067 !Local variables------------------------------- 1068 !scalars 1069 integer :: iatom,ii,iimage,imass,nmass,natom_mass 1070 real(dp) :: mtot 1071 character(len=500) :: msg 1072 !arrays 1073 real(sp) :: alea_sp(3,natom,trotter) 1074 real(dp) :: total(3) 1075 1076 !************************************************************************ 1077 1078 if (qtb_file_unit<0) then 1079 msg='QTB forces file unit should be positive!' 1080 ABI_BUG(msg) 1081 end if 1082 natom_mass=size(mass,1);nmass=size(mass,2) 1083 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1084 msg='Wrong dimensions for array mass !' 1085 ABI_BUG(msg) 1086 end if 1087 1088 !Read QTB random forces 1089 read(qtb_file_unit) alea_sp(1:3,1:natom,1:trotter) 1090 alea(:,:,:)=dble(alea_sp(:,:,:)) 1091 1092 !Make sure that the sum of random forces is zero 1093 if(zeroforce==1)then 1094 total=zero 1095 mtot=zero 1096 do iimage=1,trotter 1097 imass=min(nmass,iimage) 1098 do iatom=1,natom 1099 mtot=mtot+mass(iatom,imass) 1100 end do 1101 end do 1102 do iimage=1,trotter 1103 imass=min(nmass,iimage) 1104 do iatom=1,natom 1105 do ii=1,3 1106 total(ii)=total(ii)+langev(iatom,imass)*alea(ii,iatom,iimage) 1107 end do 1108 end do 1109 end do 1110 do iimage=1,trotter 1111 imass=min(nmass,iimage) 1112 do iatom=1,natom 1113 do ii=1,3 1114 alea(ii,iatom,iimage)= alea(ii,iatom,iimage)- & 1115 & (total(ii)*mass(iatom,imass))/(langev(iatom,imass)*mtot) 1116 end do 1117 end do 1118 end do 1119 ! now random forces have been rescaled so that their sum is zero 1120 end if 1121 1122 end subroutine pimd_langevin_random_qtb
m_pimd/pimd_mass_spring [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_mass_spring
FUNCTION
Compute masses and spring constants for PIMD. Eventually apply a coordinate transformation. Possible choices for the transformation: 0: no transformation 1: normal mode tranformation 2: staging transformation
INPUTS
inertmass(natom)=fictitious masses of atoms kt=kT constant natom=number of atoms quantummass(natom)=true masses of atoms transform=coordinate transformation: 0: no tranformation 1: normal mode transformation 2: staging transformation trotter=Trotter number
OUTPUT
SIDE EFFECTS
mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) spring(natom,mass_dim)=spring constants of atoms (mass_dim=1 or trotter)
NOTES
Back transformation (ioption=-1) not implemented !
SOURCE
2493 subroutine pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,transform,trotter) 2494 2495 !Arguments ------------------------------------ 2496 !scalars 2497 integer,intent(in) :: natom,transform,trotter 2498 real(dp),intent(in) :: kt 2499 !arrays 2500 real(dp),intent(in) :: inertmass(natom),quantummass(natom) 2501 real(dp),intent(out) :: mass(:,:),spring(:,:) 2502 !Local variables------------------------------- 2503 !scalars 2504 integer :: iimage,kk,natom_mass,natom_spring,nmass,nspring 2505 real(dp) :: gammasquare 2506 character(len=500) :: msg 2507 !arrays 2508 real(dp),allocatable :: mass_temp(:,:),lambda(:) 2509 2510 !************************************************************************ 2511 2512 natom_mass =size(mass ,1);nmass =size(mass ,2) 2513 natom_spring=size(spring,1);nspring=size(spring,2) 2514 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 2515 msg='Wrong dimensions for array mass !' 2516 ABI_BUG(msg) 2517 end if 2518 if (natom/=natom_spring.or.(nspring/=1.and.nspring/=trotter)) then 2519 msg='Wrong dimensions for array spring !' 2520 ABI_BUG(msg) 2521 end if 2522 2523 !=== No transformation =================================================== 2524 if (transform==0) then 2525 !2nd dim of mass and spring = 1 2526 mass(1:natom,1)=inertmass(1:natom) 2527 spring(1:natom,1)=quantummass(1:natom)*dble(trotter)*kt*kt 2528 2529 !=== Normal mode transformation ========================================== 2530 else if (transform==1) then 2531 2532 ABI_MALLOC(lambda,(trotter)) 2533 lambda(1)=zero; lambda(trotter)=four*dble(trotter) 2534 do kk=2,trotter/2 2535 lambda(2*kk-2)=four*dble(trotter)* & 2536 & (one-cos(two_pi*dble(kk-1)/dble(trotter))) 2537 lambda(2*kk-1)=four*dble(trotter)* & 2538 & (one-cos(two_pi*dble(kk-1)/dble(trotter))) 2539 end do 2540 2541 !normal mode masses 2542 do iimage=1,trotter 2543 mass(:,iimage)=quantummass(:)*lambda(iimage) 2544 end do 2545 2546 do iimage=1,trotter 2547 spring(:,iimage)=mass(:,iimage)*dble(trotter)*kt*kt 2548 end do 2549 2550 !fictitious masses 2551 mass(:,1)=inertmass(:) 2552 2553 !from 2 to P not changed except if adiabatic PIMD 2554 !see Hone et al, JCP 124, 154103 (2006) 2555 gammasquare=one !adiabaticity parameter 2556 do iimage=2,trotter 2557 mass(:,iimage)=mass(:,iimage)/gammasquare 2558 end do 2559 2560 ABI_FREE(lambda) 2561 2562 !=== Staging transformation ============================================== 2563 else if (transform==2) then 2564 2565 !Fictitious masses 2566 mass(1:natom,1)=inertmass(1:natom) 2567 if (nmass>1) then 2568 do iimage=2,trotter 2569 mass(1:natom,iimage)=inertmass(1:natom)*dble(iimage)/dble(iimage-1) 2570 end do 2571 end if 2572 2573 !Staging masses (mass_temp) 2574 ABI_MALLOC(mass_temp,(natom,trotter)) 2575 mass_temp(1:natom,1)=zero 2576 if (nmass>1) then 2577 do iimage=2,trotter 2578 mass_temp(1:natom,iimage)=quantummass(1:natom)*dble(iimage)/dble(iimage-1) 2579 end do 2580 end if 2581 2582 spring(1:natom,1)=mass_temp(1:natom,1)*dble(trotter)*kt*kt 2583 if (nspring>1) then 2584 do iimage=2,trotter 2585 spring(1:natom,iimage)=mass_temp(1:natom,iimage)*dble(trotter)*kt*kt 2586 end do 2587 end if 2588 ABI_FREE(mass_temp) 2589 2590 end if 2591 2592 end subroutine pimd_mass_spring
m_pimd/pimd_noseehoover_forces [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_nosehoover_forces
FUNCTION
Compute Nose-Hoover contribution to PIMD forces by adding friction force of thermostat number one
INPUTS
dzeta(3,natom,trotter,nnos)=variables of thermostats, in (atomic time unit)^(-1) used only when a coordinate transformation is applied (transfom/=0) forces(3,natom,trotter)=forces without Nose-Hoover contribution mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms nnos=number of thermostats trotter=Trotter number vel(3,natom,trotter)=velocities of atoms in each cell
OUTPUT
forces_nosehoover(3,natom,trotter)=forces including thermostat contribution
SIDE EFFECTS
SOURCE
1497 subroutine pimd_nosehoover_forces(dzeta,forces,forces_nosehoover,mass,natom,& 1498 & nnos,trotter,vel) 1499 1500 !Arguments ------------------------------------ 1501 !scalars 1502 integer,intent(in) :: natom,nnos,trotter 1503 !arrays 1504 real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),forces(3,natom,trotter) 1505 real(dp),intent(in) :: vel(3,natom,trotter) 1506 real(dp),intent(in) :: mass(:,:) 1507 real(dp),intent(out) :: forces_nosehoover(3,natom,trotter) 1508 !Local variables------------------------------- 1509 !scalars 1510 integer :: iatom,ii,iimage,imass,natom_mass,nmass 1511 character(len=500) :: msg 1512 !arrays 1513 1514 !************************************************************************ 1515 1516 natom_mass=size(mass,1);nmass=size(mass,2) 1517 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1518 msg='Wrong dimensions for array mass !' 1519 ABI_BUG(msg) 1520 end if 1521 1522 do iimage=1,trotter 1523 imass=min(nmass,iimage) 1524 do iatom=1,natom 1525 do ii=1,3 1526 forces_nosehoover(ii,iatom,iimage)=forces(ii,iatom,iimage) & 1527 & - mass(iatom,imass)*dzeta(ii,iatom,iimage,1)*vel(ii,iatom,iimage) 1528 end do 1529 end do 1530 end do 1531 1532 end subroutine pimd_nosehoover_forces
m_pimd/pimd_nosehoover_propagate [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_nosehoover_propagate
FUNCTION
Propagate thermostat variables (Nose-Hoover algorithm) - for PIMD
INPUTS
dtion=time step dzeta(3,natom,trotter,nnos)=time derivative of zeta at time t itimimage=index of time step mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms nnos=number of thermostats qmass(nnos)=masses of thermostats temperature=temperature transform=coordinate transformation: 0: no tranformation 1: normal mode transformation 2: staging transformation trotter=Trotter number vel(3,natom,trotter)=velocities of atoms in each cell zeta(3,natom,trotter,nnos)=variables of thermostats, in (atomic time unit)^(-1) used only when no coordinate transformation is applied (transfom==0) at time t zeta_prev(3,natom,trotter,nnos)=previous value of zeta (t-dt)
OUTPUT
zeta_next(3,natom,trotter,nnos)=next value of zeta (t+dt)
SIDE EFFECTS
SOURCE
1896 subroutine pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,temperature,& 1897 & trotter,vel,zeta,zeta_next,zeta_prev,itimimage,transform) 1898 1899 !Arguments ------------------------------------ 1900 !scalars 1901 integer,intent(in) :: itimimage,natom,nnos,transform,trotter 1902 real(dp),intent(in) :: dtion,temperature 1903 !arrays 1904 real(dp),intent(in) :: qmass(nnos),vel(3,natom,trotter) 1905 real(dp),intent(in) :: mass(:,:) 1906 real(dp),intent(in) :: dzeta(3,natom,trotter,nnos),zeta_prev(3,natom,trotter,nnos) 1907 real(dp),intent(in) :: zeta(3,natom,trotter,nnos) 1908 real(dp),intent(out) :: zeta_next(3,natom,trotter,nnos) 1909 !Local variables------------------------------- 1910 !scalars 1911 integer :: iatom,ii,iimage,inos,natom_mass,nmass 1912 real(dp) :: kt 1913 character(len=500) :: msg 1914 !arrays 1915 real(dp),allocatable :: thermforces(:,:,:,:) 1916 1917 !************************************************************************ 1918 1919 natom_mass=size(mass,1);nmass=size(mass,2) 1920 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1921 msg='Wrong dimensions for array mass !' 1922 ABI_BUG(msg) 1923 end if 1924 if (nnos<3) then 1925 msg='Not available for nnos<3 !' 1926 ABI_BUG(msg) 1927 end if 1928 1929 kt=temperature*kb_HaK 1930 ABI_MALLOC(thermforces,(3,natom,trotter,nnos)) 1931 1932 !forces on the thermostats 1933 do inos=1,nnos 1934 if(inos==1)then 1935 do iimage=1,trotter 1936 do iatom=1,natom 1937 do ii=1,3 1938 select case(transform) 1939 case(0) 1940 thermforces(ii,iatom,iimage,1)=mass(iatom,1)*(vel(ii,iatom,iimage)**2)-kt 1941 case(1) 1942 thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt 1943 case(2) 1944 thermforces(ii,iatom,iimage,1)=mass(iatom,iimage)*(vel(ii,iatom,iimage)**2)-kt 1945 end select 1946 end do 1947 end do 1948 end do 1949 elseif(inos==(nnos-1))then 1950 do iimage=1,trotter 1951 do iatom=1,natom 1952 do ii=1,3 1953 thermforces(ii,iatom,iimage,inos)=& 1954 & (qmass(nnos-2)*(dzeta(ii,iatom,iimage,nnos-2)**2)-kt) + & 1955 & (qmass(nnos )*(dzeta(ii,iatom,iimage,nnos) **2)-kt) 1956 end do 1957 end do 1958 end do 1959 elseif(inos==nnos)then 1960 do iimage=1,trotter 1961 do iatom=1,natom 1962 do ii=1,3 1963 thermforces(ii,iatom,iimage,inos)=& 1964 & qmass(nnos-1)*(dzeta(ii,iatom,iimage,nnos-1)**2)-kt 1965 end do 1966 end do 1967 end do 1968 else 1969 do iimage=1,trotter 1970 do iatom=1,natom 1971 do ii=1,3 1972 thermforces(ii,iatom,iimage,inos)=& 1973 & qmass(inos-1)*(dzeta(ii,iatom,iimage,inos-1)**2)-kt 1974 end do 1975 end do 1976 end do 1977 end if 1978 end do 1979 1980 select case(itimimage) 1981 case(1) !taylor 1982 1983 do inos=1,nnos 1984 if(inos==1)then 1985 zeta_next(:,:,:,1)=zeta(:,:,:,1)+ dzeta(:,:,:,1)*dtion + & 1986 & (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* & 1987 & dtion*dtion/(two*qmass(1)) 1988 elseif(inos==(nnos-1))then 1989 zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + & 1990 & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* & 1991 & dtion*dtion/(two*qmass(inos)) 1992 elseif(inos==nnos)then 1993 zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + & 1994 & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* & 1995 & dtion*dtion/(two*qmass(inos)) 1996 else 1997 zeta_next(:,:,:,inos)=zeta(:,:,:,inos)+ dzeta(:,:,:,inos)*dtion + & 1998 & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* & 1999 & dtion*dtion/(two*qmass(inos)) 2000 end if 2001 end do 2002 2003 case default !verlet 2004 2005 do inos=1,nnos 2006 if(inos==1)then 2007 zeta_next(:,:,:,1)=two*zeta(:,:,:,1) - zeta_prev(:,:,:,1) + & 2008 & (thermforces(:,:,:,1)-qmass(1)*dzeta(:,:,:,1)*dzeta(:,:,:,2))* & 2009 & dtion*dtion/qmass(1) 2010 elseif(inos==(nnos-1))then 2011 zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + & 2012 & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,nnos))* & 2013 & dtion*dtion/qmass(inos) 2014 elseif(inos==nnos)then 2015 zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + & 2016 & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,nnos-1)*dzeta(:,:,:,nnos))* & 2017 & dtion*dtion/qmass(inos) 2018 else 2019 zeta_next(:,:,:,inos)=two*zeta(:,:,:,inos) - zeta_prev(:,:,:,inos) + & 2020 & (thermforces(:,:,:,inos)-qmass(inos)*dzeta(:,:,:,inos)*dzeta(:,:,:,inos+1))* & 2021 & dtion*dtion/qmass(inos) 2022 end if 2023 end do 2024 2025 end select 2026 2027 ABI_FREE(thermforces) 2028 2029 end subroutine pimd_nosehoover_propagate
m_pimd/pimd_nullify [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_nullify
FUNCTION
Nullify the content of a datastructure of type pimd_type.
INPUTS
OUTPUT
SIDE EFFECTS
pimd_param=datastructure of type pimd_type. several parameters for Path-Integral MD.
SOURCE
219 subroutine pimd_nullify(pimd_param) 220 221 implicit none 222 223 !Arguments ------------------------------------ 224 !scalars 225 type(pimd_type),intent(inout) :: pimd_param 226 227 !************************************************************************ 228 229 pimd_param%adpimd = 0 230 pimd_param%constraint = 0 231 pimd_param%irandom = -1 232 pimd_param%nnos = -1 233 pimd_param%ntypat = -1 234 pimd_param%optcell = -1 235 pimd_param%pitransform = -1 236 pimd_param%qtb_file_unit= -1 237 pimd_param%traj_unit = -1 238 pimd_param%use_qtb = 0 239 pimd_param%adpimd_gamma = one 240 pimd_param%vis = zero 241 pimd_param%bmass = zero 242 pimd_param%dtion = zero 243 pimd_param%friction = zero 244 nullify(pimd_param%mdtemp) 245 nullify(pimd_param%pimass) 246 nullify(pimd_param%strtarget) 247 nullify(pimd_param%amu) 248 nullify(pimd_param%qmass) 249 nullify(pimd_param%typat) 250 nullify(pimd_param%wtatcon) 251 252 end subroutine pimd_nullify
m_pimd/pimd_predict_taylor [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_predict_taylor
FUNCTION
Predict new atomic positions using a Taylor algorithm (first time step) - for PIMD
INPUTS
dtion=time step forces(3,natom,trotter)=PIMD forces on atoms in each cell mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms trotter=Trotter number vel(3,natom,trotter)=velocities of atoms in each cell xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t
OUTPUT
xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt
SIDE EFFECTS
SOURCE
1757 subroutine pimd_predict_taylor(dtion,forces,mass,natom,trotter,vel,xcart,xcart_next) 1758 1759 !Arguments ------------------------------------ 1760 !scalars 1761 integer,intent(in) :: natom,trotter 1762 real(dp),intent(in) :: dtion 1763 !arrays 1764 real(dp),intent(in) :: forces(3,natom,trotter),vel(3,natom,trotter),xcart(3,natom,trotter) 1765 real(dp),intent(in) :: mass(:,:) 1766 real(dp),intent(out) :: xcart_next(3,natom,trotter) 1767 !Local variables------------------------------- 1768 !scalars 1769 integer :: iatom,ii,iimage,imass,natom_mass,nmass 1770 character(len=500) :: msg 1771 !arrays 1772 1773 !************************************************************************ 1774 1775 natom_mass=size(mass,1);nmass=size(mass,2) 1776 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1777 msg='Wrong dimensions for array mass !' 1778 ABI_BUG(msg) 1779 end if 1780 1781 do iimage=1,trotter 1782 imass=min(nmass,iimage) 1783 do iatom=1,natom 1784 do ii=1,3 1785 xcart_next(ii,iatom,iimage)=xcart(ii,iatom,iimage) & 1786 & + half*dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass) & 1787 & + dtion*vel(ii,iatom,iimage) 1788 end do 1789 end do 1790 end do 1791 1792 end subroutine pimd_predict_taylor
m_pimd/pimd_predict_verlet [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_predict_verlet
FUNCTION
Predict new atomic positions using a Verlet algorithm - for PIMD
INPUTS
dtion=time step forces(3,natom,trotter)=PIMD forces on atoms in each cell mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms trotter=Trotter number xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t xcart_prev(3,natom,trotter)=cartesian coordinates of atoms in each cell at t-dt
OUTPUT
xcart_next(3,natom,trotter)=cartesian coordinates of atoms in each cell at t+dt
SIDE EFFECTS
SOURCE
1820 subroutine pimd_predict_verlet(dtion,forces,mass,natom,trotter,xcart,xcart_next,xcart_prev) 1821 1822 !Arguments ------------------------------------ 1823 !scalars 1824 integer,intent(in) :: natom,trotter 1825 real(dp),intent(in) :: dtion 1826 !arrays 1827 real(dp),intent(in) :: forces(3,natom,trotter) 1828 real(dp),intent(in) :: xcart(3,natom,trotter),xcart_prev(3,natom,trotter) 1829 real(dp),intent(in) :: mass(:,:) 1830 real(dp),intent(out) :: xcart_next(3,natom,trotter) 1831 !Local variables------------------------------- 1832 !scalars 1833 integer :: iatom,ii,iimage,imass,natom_mass,nmass 1834 character(len=500) :: msg 1835 !arrays 1836 1837 !************************************************************************ 1838 1839 natom_mass=size(mass,1);nmass=size(mass,2) 1840 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1841 msg='Wrong dimensions for array mass !' 1842 ABI_BUG(msg) 1843 end if 1844 1845 do iimage=1,trotter 1846 imass=min(nmass,iimage) 1847 do iatom=1,natom 1848 do ii=1,3 1849 xcart_next(ii,iatom,iimage)= & 1850 & two*xcart(ii,iatom,iimage) & 1851 & - xcart_prev(ii,iatom,iimage) & 1852 & + dtion*dtion*forces(ii,iatom,iimage)/mass(iatom,imass) 1853 end do 1854 end do 1855 end do 1856 1857 end subroutine pimd_predict_verlet
m_pimd/pimd_print [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_print
FUNCTION
Print out results related to PIMD (for given time step)
INPUTS
constraint=type of constraint eventually applied (on a reaction coordinate) constraint_output(2)=several (real) data to be output when a constraint has been applied eharm=harmonic energy eharm_virial=harmonic energy from virial epot=potential energy forces(3,natom,trotter)=forces on atoms in each cell inertmass(natom)=inertial masses of atoms irestart=1 if this is a restart itimimage=index of time step kt=Thermal energy K_b.T natom=number of atoms optcell=option for cell evolution prtstress=flag for stress tensor printing prtvolimg=printing volume rprim(3,3)=dimensionless real space primitive translations stress(3,3,3)=stress tensor (first dimension corresponds to 3 different estimators) temperature1,temperature2=temperatures at t and t+dt traj_unit=flag activating printing of the trajectory in an external file\ if >0, indicates the unit number of the trajectory file trotter=Trotter number vel(3,natom,trotter)=velocities of atoms in in each cell vel_cell(3,3)= velocities of cell parameters (time derivative of rprimd) xcart(3,natom,trotter)=cartesian coordinates of atoms in in each cell xred(3,natom,trotter)=reduced coordinates of atoms in in each cell
OUTPUT
-- only printing --
SIDE EFFECTS
SOURCE
572 subroutine pimd_print(constraint,constraint_output,eharm,eharm_virial,epot,& 573 & forces,inertmass,irestart,itimimage,kt,natom,optcell,prtstress,& 574 & prtvolimg,rprimd,stress,temperature1,temperature2,traj_unit,& 575 & trotter,vel,vel_cell,xcart,xred) 576 577 !Arguments ------------------------------------ 578 !scalars 579 integer,intent(in) :: constraint,irestart,itimimage,natom,optcell 580 integer,intent(in) :: prtstress,prtvolimg,traj_unit,trotter 581 real(dp),intent(in) :: eharm,eharm_virial,epot,kt,temperature1,temperature2 582 !arrays 583 real(dp),intent(in) :: constraint_output(2) 584 real(dp),intent(in) :: forces(3,natom,trotter),inertmass(natom) 585 real(dp),intent(in) :: rprimd(3,3),stress(3,3,3),vel(3,natom,trotter),vel_cell(3,3) 586 real(dp),intent(in) :: xcart(3,natom,trotter),xred(3,natom,trotter) 587 !Local variables------------------------------- 588 !scalars 589 integer :: iatom,ii,iimage 590 real(dp) :: mtot 591 character(len=500) :: msg 592 !arrays 593 real(dp) :: acell(3),cdm_cart(3),cdm_red(3),forcetot(3),rprim(3,3) 594 real(dp),allocatable :: centroid(:,:),qudeloc(:) 595 596 !************************************************************************ 597 598 !Temperature 599 if(itimimage==1)then 600 msg=ch10 601 if(mod(irestart,10)==0) then 602 write(msg,'(2a)') ch10,' This is a PIMD calculation from scratch' 603 else if (mod(irestart,10)==1) then 604 write(msg,'(2a)') ch10,' This is a RESTART calculation' 605 end if 606 call wrtout(ab_out,msg,'COLL') 607 call wrtout(std_out,msg,'COLL') 608 write(msg,'(a,f12.5,a)') & 609 & ' In the initial configuration, the temperature is ',temperature1,' K' 610 call wrtout(ab_out,msg,'COLL') 611 call wrtout(std_out,msg,'COLL') 612 end if 613 write(msg,'(2a,i5,a,f12.5,a)') ch10,& 614 & ' At PIMD time step ',itimimage,', the temperature is',temperature2,' K' 615 call wrtout(ab_out,msg,'COLL') 616 call wrtout(std_out,msg,'COLL') 617 618 !Energies 619 write(msg,'(4a,f18.9,a,a,a,f18.9,a,a)') ch10,& 620 & ' Energy:',ch10, & 621 & ' Internal energy (PRIMITIVE estimator) =',onehalf*dble(natom*trotter)*kt-eharm+epot ,' Ha',ch10, & 622 & ' Internal energy (VIRIAL estimator) =',onehalf*dble(natom)*kt+eharm_virial+epot,' Ha',ch10 623 call wrtout(ab_out,msg,'COLL') 624 call wrtout(std_out,msg,'COLL') 625 626 !Stress tensor and pressure 627 write(msg,'(2a,3(2a,3f18.9))') ch10,& 628 & ' Stress tensor from PRIMITIVE estimator (Ha/Bohr^3):',ch10, & 629 & ' ',stress(1,1,1),stress(1,1,2),stress(1,1,3),ch10, & 630 & ' ',stress(1,2,1),stress(1,2,2),stress(1,2,3),ch10, & 631 & ' ',stress(1,3,1),stress(1,3,2),stress(1,3,3) 632 if (prtstress==1) then 633 call wrtout(ab_out,msg,'COLL') 634 end if 635 call wrtout(std_out,msg,'COLL') 636 write(msg,'(a,f18.9,a)') ' Pressure (primitive estimator) =', & 637 & -third*(stress(1,1,1)+stress(1,2,2)+stress(1,3,3))*HaBohr3_GPa,' GPa' 638 call wrtout(ab_out,msg,'COLL') 639 call wrtout(std_out,msg,'COLL') 640 641 !Data related to constraint eventually applied 642 if (constraint/=0) then 643 if (constraint==1) write(msg,'(2a)') ch10,' Blue Moon Ensemble method is activated:' 644 call wrtout(ab_out,msg,'COLL') 645 call wrtout(std_out,msg,'COLL') 646 write(msg,'(a,f18.10,2a,f18.10)') & 647 & ' - Reaction coordinate =',constraint_output(1),ch10,& 648 & ' - Instantaneous force on the reaction coord. =',constraint_output(2) 649 call wrtout(ab_out,msg,'COLL') 650 call wrtout(std_out,msg,'COLL') 651 end if 652 653 !Total force 654 if (prtvolimg<=1) then 655 forcetot=zero 656 do iimage=1,trotter 657 do iatom=1,natom 658 do ii=1,3 659 forcetot(ii)=forcetot(ii)+forces(ii,iatom,iimage) 660 end do 661 end do 662 end do 663 write(msg,'(2a,3f18.10)') ch10,' Total force=',forcetot(1:3) 664 call wrtout(std_out,msg,'COLL') 665 end if 666 667 !position of mass center 668 if (prtvolimg<=1) then 669 mtot=zero;cdm_cart=zero;cdm_red=zero 670 do iimage=1,trotter 671 do iatom=1,natom 672 cdm_cart(:)=cdm_cart(:)+inertmass(iatom)*xcart(:,iatom,iimage) 673 cdm_red (:)=cdm_red (:)+inertmass(iatom)*xred (:,iatom,iimage) 674 mtot=mtot+inertmass(iatom) 675 end do 676 end do 677 cdm_cart=cdm_cart/mtot 678 cdm_red =cdm_red /mtot 679 write(msg,'(3a,3x,3f18.10,3a,3x,3f18.10)') ch10,& 680 & ' Center of mass, in cartes. coordinates :',ch10,cdm_cart(:),ch10,& 681 & ' Center of mass, in reduced coordinates :',ch10,cdm_red(:) 682 call wrtout(std_out,msg,'COLL') 683 call wrtout(ab_out,msg,'COLL') 684 end if 685 686 !Positions 687 write(msg,'(2a)') ch10,' Atomic positions:' 688 call wrtout(std_out,msg,'COLL') 689 call wrtout(ab_out,msg,'COLL') 690 do iimage=1,trotter 691 select case(iimage) 692 case(1) 693 write(msg,'(a)') ' xred' 694 case(2,3,4,5,6,7,8,9) 695 write(msg,'(a,i1,a)') ' xred_',iimage,'img' 696 case(10:99) 697 write(msg,'(a,i2,a)') ' xred_',iimage,'img' 698 case default 699 write(msg,'(a,i3,a)') ' xred_',iimage,'img' 700 end select 701 call wrtout(std_out,msg,'COLL') 702 call wrtout(ab_out,msg,'COLL') 703 if (traj_unit>0) then 704 call wrtout(traj_unit,msg,'COLL') 705 end if 706 do iatom=1,natom 707 write(msg,'(3f18.10)') xred(1:3,iatom,iimage) 708 call wrtout(std_out,msg,'COLL') 709 call wrtout(ab_out,msg,'COLL') 710 if (traj_unit>0) then 711 call wrtout(traj_unit,msg,'COLL') 712 end if 713 end do 714 end do 715 716 !Velocities 717 write(msg,'(2a)') ch10,' Velocities:' 718 call wrtout(std_out,msg,'COLL') 719 call wrtout(ab_out,msg,'COLL') 720 do iimage=1,trotter 721 select case(iimage) 722 case(1) 723 write(msg,'(a)') ' vel' 724 case(2,3,4,5,6,7,8,9) 725 write(msg,'(a,i1,a)') ' vel_',iimage,'img' 726 case(10:99) 727 write(msg,'(a,i2,a)') ' vel_',iimage,'img' 728 case default 729 write(msg,'(a,i3,a)') ' vel_',iimage,'img' 730 end select 731 call wrtout(std_out,msg,'COLL') 732 call wrtout(ab_out,msg,'COLL') 733 if (traj_unit>0) then 734 call wrtout(traj_unit,msg,'COLL') 735 end if 736 do iatom=1,natom 737 write(msg,'(3f18.10)') vel(1:3,iatom,iimage) 738 call wrtout(std_out,msg,'COLL') 739 call wrtout(ab_out,msg,'COLL') 740 if (traj_unit>0) then 741 call wrtout(traj_unit,msg,'COLL') 742 end if 743 end do 744 end do 745 746 if (optcell>0) then 747 748 call mkradim(acell,rprim,rprimd) 749 750 ! Time derivative of rprimd 751 write(msg,'(2a)') ch10,' Time derivative of rprimd:' 752 call wrtout(std_out,msg,'COLL') 753 call wrtout(ab_out,msg,'COLL') 754 write(msg,'(2a,3(3f18.10))') ' vel_cell',ch10,vel_cell(:,:) 755 call wrtout(std_out,msg,'COLL') 756 call wrtout(ab_out,msg,'COLL') 757 if (traj_unit>0) then 758 call wrtout(traj_unit,msg,'COLL') 759 end if 760 761 ! rprimd 762 write(msg,'(2a)') ch10,' Cell parameters:' 763 call wrtout(std_out,msg,'COLL') 764 call wrtout(ab_out,msg,'COLL') 765 write(msg,'(2a,3(3f18.10),3a,3f18.10)') ' rprim',ch10,rprim(:,:),ch10,& 766 & ' acell',ch10,acell(:) 767 call wrtout(std_out,msg,'COLL') 768 call wrtout(ab_out,msg,'COLL') 769 if (traj_unit>0) then 770 call wrtout(traj_unit,msg,'COLL') 771 end if 772 773 end if 774 775 !Centroids and wave-packet spatial spreads 776 ABI_MALLOC(centroid,(3,natom)) 777 ABI_MALLOC(qudeloc,(natom)) 778 centroid=zero;qudeloc=zero 779 do iimage=1,trotter 780 do iatom=1,natom 781 do ii=1,3 782 centroid(ii,iatom)=centroid(ii,iatom)+xcart(ii,iatom,iimage) 783 end do 784 end do 785 end do 786 centroid=centroid/dble(trotter) 787 do iimage=1,trotter 788 do iatom=1,natom 789 do ii=1,3 790 qudeloc(iatom)=qudeloc(iatom)+((xcart(ii,iatom,iimage)-centroid(ii,iatom))**2) 791 end do 792 end do 793 end do 794 qudeloc(:)=sqrt(qudeloc(:)/dble(trotter)) 795 write(msg,'(4a)') ch10,' Centroids and wave-packet spatial spreads (cart. coord.):',ch10,& 796 & ' iat centroid_x centroid_y centroid_z spatial_spread' 797 call wrtout(std_out,msg,'COLL') 798 do iatom=1,natom 799 write(msg,'(i4,4f18.10)') iatom,centroid(1:3,iatom),qudeloc(iatom) 800 call wrtout(std_out,msg,'COLL') 801 end do 802 ABI_FREE(centroid) 803 ABI_FREE(qudeloc) 804 805 !Fake statement 806 return;ii=prtvolimg 807 808 end subroutine pimd_print
m_pimd/pimd_skip_qtb [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_skip_qtb
FUNCTION
Only relevant in case of PI-QTB: Skip a line in a QTB random force file
INPUTS
pimd_param=datastructure of type pimd_type. several parameters for Path-Integral MD.
OUTPUT
SOURCE
398 subroutine pimd_skip_qtb(pimd_param) 399 400 !Arguments ------------------------------------ 401 !scalars 402 type(pimd_type),intent(in) :: pimd_param 403 !arrays 404 !Local variables------------------------------- 405 !scalars 406 character(len=500) :: msg 407 !arrays 408 409 !************************************************************************ 410 411 if (pimd_param%use_qtb==0) return 412 413 if (pimd_param%qtb_file_unit<0) then 414 msg='QTB forces file unit should be positive!' 415 ABI_BUG(msg) 416 end if 417 418 !Skip one line QTB random forces file 419 read(pimd_param%qtb_file_unit) 420 421 end subroutine pimd_skip_qtb
m_pimd/pimd_stresses [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_stresses
FUNCTION
In the case od PIMD, compute the pressure tensor from virial theorem
INPUTS
mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) natom=number of atoms quantummass(natom)=quantum masses of atoms stressin(3,3,trotter)=electronic stress tensor for each image temperature=temperature (could be instantaneous temp. or thermostat temp.) temperature_therm=thermostat temperature trotter=Trotter number vel(3,natom,trotter)=velocities of atoms in each cell volume=volume of each cell (common to all cells) xcart(3,natom,trotter)=cartesian coordinates of atoms in each cell at t temperature=thermostat temperature
OUTPUT
stress_pimd(3,3,3)=stress tensor for PIMD First dimension (3) corresponds to 3 different pressure estimators
SIDE EFFECTS
SOURCE
1565 subroutine pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,& 1566 & temperature,temperature_therm,trotter,vel,volume,xcart) 1567 1568 !Arguments ------------------------------------ 1569 !scalars 1570 integer,intent(in) :: natom,trotter 1571 real(dp),intent(in) :: temperature,temperature_therm,volume 1572 !arrays 1573 real(dp),intent(in) :: vel(3,natom,trotter),xcart(3,natom,trotter) 1574 real(dp),intent(in) :: mass(:,:),stressin(3,3,trotter),quantummass(natom) 1575 real(dp),intent(out) :: stress_pimd(3,3,3) 1576 !Local variables------------------------------- 1577 !scalars 1578 integer :: iatom,ii,iimage,imass,jj,natom_mass,nmass,iimagep 1579 real(dp) :: stress_tmp(3,3,3),omega2,kt,kt_therm 1580 character(len=500) :: msg 1581 1582 !************************************************************************ 1583 1584 natom_mass=size(mass,1);nmass=size(mass,2) 1585 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=trotter)) then 1586 msg='Wrong dimensions for array mass !' 1587 ABI_BUG(msg) 1588 end if 1589 1590 kt=temperature*kb_HaK 1591 kt_therm=temperature_therm*kb_HaK 1592 stress_pimd=zero 1593 1594 !I-PRIMITIVE ESTIMATOR 1595 !1-Kinetic part 1596 do ii=1,3 1597 stress_pimd(1,ii,ii)=dble(natom)*dble(trotter)*kt/volume 1598 end do 1599 1600 !2-Potential part 1601 do iimage=1,trotter 1602 do ii=1,3 1603 do jj=1,3 1604 stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-stressin(ii,jj,iimage)/dble(trotter) 1605 !minus to convert stress into pressure 1606 end do 1607 end do 1608 end do 1609 1610 !3-Contribution from springs 1611 omega2=dble(trotter)*kt_therm*kt_therm 1612 do iimage=1,trotter 1613 iimagep=iimage+1 1614 if(iimage==trotter) then 1615 iimagep=1 1616 end if 1617 do iatom=1,natom 1618 do ii=1,3 1619 do jj=1,3 1620 stress_pimd(1,ii,jj)=stress_pimd(1,ii,jj)-quantummass(iatom)*omega2* & 1621 & (xcart(ii,iatom,iimagep)-xcart(ii,iatom,iimage))* & 1622 & (xcart(jj,iatom,iimagep)-xcart(jj,iatom,iimage))/volume 1623 end do 1624 end do 1625 end do 1626 end do 1627 1628 !II-Average of classical pressures 1629 !1-Kinetic part 1630 do iimage=1,trotter 1631 imass=min(nmass,iimage) 1632 do iatom=1,natom 1633 do ii=1,3 1634 do jj=1,3 1635 stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)+ & 1636 & mass(iatom,imass)*vel(ii,iatom,iimage)*vel(jj,iatom,iimage)/volume 1637 end do 1638 end do 1639 end do 1640 end do 1641 1642 !2-Contribution from electronic stress 1643 do iimage=1,trotter 1644 do ii=1,3 1645 do jj=1,3 1646 stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)-stressin(ii,jj,iimage) 1647 end do 1648 end do 1649 end do 1650 do ii=1,3 1651 do jj=1,3 1652 stress_pimd(2,ii,jj)=stress_pimd(2,ii,jj)/dble(trotter) 1653 end do 1654 end do 1655 1656 !III-pressure from VIRIAL estimator: stress_pimd(3,:,:) 1657 !1-kinetic part 1658 ! do ii=1,3 1659 ! stress_pimd(3,ii,ii)=dfloat(natom)*kt/volume 1660 ! end do 1661 1662 !Symmetrize internal pressure 1663 stress_tmp=stress_pimd 1664 stress_pimd(:,2,1)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1)) 1665 stress_pimd(:,1,2)=half*(stress_tmp(:,1,2)+stress_tmp(:,2,1)) 1666 stress_pimd(:,3,1)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1)) 1667 stress_pimd(:,1,3)=half*(stress_tmp(:,1,3)+stress_tmp(:,3,1)) 1668 stress_pimd(:,2,3)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3)) 1669 stress_pimd(:,3,2)=half*(stress_tmp(:,3,2)+stress_tmp(:,2,3)) 1670 1671 end subroutine pimd_stresses
m_pimd/pimd_temperature [ Functions ]
[ Top ] [ m_pimd ] [ Functions ]
NAME
pimd_temperature
FUNCTION
Compute temperature from velocities and masses
INPUTS
mass(natom,mass_dim)=masses of atoms (mass_dim=1 or trotter) vel(3,natom,nimage)=velocities for each image of the cell
OUTPUT
pimd_temperature=temperature (from all images of the cell)
SIDE EFFECTS
SOURCE
489 function pimd_temperature(mass,vel) 490 491 !Arguments ------------------------------------ 492 !scalars 493 real(dp) :: pimd_temperature 494 !arrays 495 real(dp),intent(in) :: mass(:,:),vel(:,:,:) 496 !Local variables------------------------------- 497 !scalars 498 integer :: iatom,idir,iimage,imass,natom,natom_mass,ndir,nimage,nmass 499 real(dp) :: v2 500 character(len=500) :: msg 501 !arrays 502 503 !************************************************************************ 504 505 ndir=size(vel,1);natom=size(vel,2);nimage=size(vel,3) 506 natom_mass=size(mass,1);nmass=size(mass,2) 507 if (ndir/=3.or.natom<=0.or.nimage<=0) then 508 msg='Wrong sizes for vel array !' 509 ABI_BUG(msg) 510 end if 511 if (natom/=natom_mass.or.(nmass/=1.and.nmass/=nimage)) then 512 msg='Wrong dimensions for array mass !' 513 ABI_BUG(msg) 514 end if 515 516 v2=zero 517 do iimage=1,nimage 518 imass=min(nmass,iimage) 519 do iatom=1,natom 520 do idir=1,3 521 v2=v2+vel(idir,iatom,iimage)*vel(idir,iatom,iimage)*mass(iatom,imass) 522 end do 523 end do 524 end do 525 pimd_temperature=v2/(dble(3*natom*nimage)*kb_HaK) 526 527 end function pimd_temperature
m_pimd/pimd_type [ Types ]
NAME
pimd_type
FUNCTION
Datatype with the variables required to perform PIMD
NOTES
SOURCE
80 type,public :: pimd_type 81 ! Scalars 82 integer :: adpimd 83 integer :: constraint 84 integer :: irandom 85 integer :: nnos 86 integer :: ntypat 87 integer :: optcell 88 integer :: pitransform 89 integer :: traj_unit 90 integer :: use_qtb 91 integer :: qtb_file_unit 92 real(dp) :: adpimd_gamma 93 real(dp) :: vis 94 real(dp) :: bmass 95 real(dp) :: dtion 96 real(dp) :: friction 97 ! Arrays 98 integer ,pointer :: typat(:) ! This pointer is associated with dtset%typat 99 real(dp),pointer :: amu(:) ! This pointer is associated with dtset%%amu_orig(:,1) 100 real(dp),pointer :: mdtemp(:) ! This pointer is associated with dtset%mdtemp 101 real(dp),pointer :: pimass(:) ! This pointer is associated with dtset%pimass 102 real(dp),pointer :: qmass(:) ! This pointer is associated with dtset%qmass 103 real(dp),pointer :: strtarget(:) ! This pointer is associated with dtset%strtarget 104 real(dp),pointer :: wtatcon(:,:,:) ! This pointer is associated with dtset%wtatcon 105 real(dp),allocatable :: zeta_prev(:,:,:,:) 106 real(dp),allocatable :: zeta (:,:,:,:) 107 real(dp),allocatable :: zeta_next(:,:,:,:) 108 real(dp),allocatable :: dzeta (:,:,:,:) 109 end type pimd_type