Skip to content

Commit

Permalink
PIMD implemented properly & tests added
Browse files Browse the repository at this point in the history
  • Loading branch information
@mbarnfield63
mbarnfield63 committed Mar 11, 2024
1 parent e2a8d7b commit a8da962
Show file tree
Hide file tree
Showing 8 changed files with 1,596 additions and 8 deletions.
14 changes: 6 additions & 8 deletions src/force_h2o.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads)
integer, intent(in) :: natom, nbeads
! Internal water coordinates
real(DP) :: rOH1, rOH2, aHOH_rad
real(DP) :: rij(1, 3)
real(DP) :: rij(nbeads, 3)
real(DP) :: Epot(nbeads)
integer :: iw

Expand Down Expand Up @@ -125,7 +125,7 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
! Schwenke calculated peterbed geometry energy
real(DP) :: Epot_delta(1)

real(DP) :: Eclas_orig, Eclas_plus
real(DP) :: Eclas_orig
real(DP) :: delta = 5.0E-5_DP
integer :: i, j, k

Expand Down Expand Up @@ -159,18 +159,16 @@ subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)

new_rij(1, :) = [new_rOH1, new_rOH2, new_aHOH_rad]

call h2o_pot_schwenke(new_rij, Epot_delta(1), nbeads)

Eclas_plus = Epot_delta(1)
call h2o_pot_schwenke(new_rij, Epot_delta(1), 1)

! Calculate the numerical force
select case (k)
case (1)
fx(i, j) = -(Eclas_plus - Eclas_orig) / delta
fx(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
case (2)
fy(i, j) = -(Eclas_plus - Eclas_orig) / delta
fy(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
case (3)
fz(i, j) = -(Eclas_plus - Eclas_orig) / delta
fz(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
end select
end do
end do
Expand Down
2 changes: 2 additions & 0 deletions tests/H2O_SCHWENKE_PIMD/est_energy.dat.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# Time[fs] E-potential E-primitive E-virial CumulAvg_prim CumulAvg_vir
0.48 0.2032534245E-03 0.4240025844E-01 0.4610886063E-02 0.0000000000E+00 0.0000000000E+00
32 changes: 32 additions & 0 deletions tests/H2O_SCHWENKE_PIMD/forces.xyz.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
3
net force: -0.34900E-04 -0.21872E-04 -0.23623E-06 torque force: -0.11636E-06 0.11789E-06 0.12688E-04
O 0.1406278399E-02 -0.7643573452E-04 -0.4539698562E-05
H -0.4915237113E-03 0.5916075151E-03 0.7610904203E-05
H -0.9183288214E-03 -0.5174496949E-03 -0.3070025751E-05
O 0.1029504537E-02 -0.3137513875E-03 -0.7261943491E-05
H -0.1096923748E-02 0.5550899357E-03 0.9241624504E-05
H 0.6388830400E-04 -0.2434761981E-03 -0.1999998591E-05
O 0.8080863472E-03 -0.1105513313E-02 -0.6980333922E-05
H -0.1174807221E-02 0.8695034946E-03 0.9433782651E-05
H 0.3633443635E-03 0.2338748312E-03 -0.2503661399E-05
O -0.1152977997E-02 0.1000189522E-03 0.6741341207E-05
H 0.4578461668E-03 -0.4920100495E-03 -0.2205022587E-05
H 0.6915218342E-03 0.3897813730E-03 -0.4534861083E-05
O -0.2491688014E-02 -0.9286233146E-03 -0.3212782834E-04
H 0.1887885553E-03 -0.4651981520E-03 -0.8959672607E-05
H 0.2299520787E-02 0.1391728496E-02 0.4103745521E-04
O -0.6996351407E-03 -0.1407817375E-02 0.2441041569E-04
H -0.4892135576E-03 0.4127463371E-03 -0.8007478791E-05
H 0.1185525436E-02 0.9929424867E-03 -0.1646446403E-04
O 0.1156527630E-03 -0.2174677061E-04 0.1258084470E-05
H -0.9827495886E-04 0.5499903011E-04 -0.1533251695E-05
H -0.2096749453E-04 -0.3547954388E-04 0.2731750706E-06
O -0.1569677876E-03 -0.6794758758E-03 0.1410986050E-04
H -0.2039187457E-03 0.2854735131E-03 -0.6374207251E-05
H 0.3574246109E-03 0.3917950652E-03 -0.7764241868E-05
O -0.4446000823E-04 0.1103340000E-03 -0.9740871998E-06
H 0.3852600226E-04 -0.7298840577E-04 0.7055280500E-06
H 0.2318156237E-05 -0.3957809077E-04 0.2724330799E-06
O -0.1330421239E-02 -0.7888721397E-03 0.7259842394E-05
H 0.4495847032E-03 -0.1294748862E-03 0.2919610472E-05
H 0.8773974682E-03 0.9161235934E-03 -0.1020951249E-04
29 changes: 29 additions & 0 deletions tests/H2O_SCHWENKE_PIMD/input.in
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
&general
nstep=1,
irest=0,
idebug=0
iknow=0

pot='_h2o_'
ipimd=1,
nwalk=10,
istage=1,
nproc=1

dt=20.,
irandom=131313,

nwrite=1,
nwritex=1,
nrest=1,
nwritef=1, ! write forces
/

&nhcopt
inose=1,
temp=298.15,
tau0=0.0015
nrespnose=3
nyosh=7
rem_comrot=.false.
/
5 changes: 5 additions & 0 deletions tests/H2O_SCHWENKE_PIMD/mini.xyz
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
3

O 0.000 0.118 0.000
H 0.757 -0.472 0.000
H -0.757 -0.472 0.000
50 changes: 50 additions & 0 deletions tests/H2O_SCHWENKE_PIMD/movie.xyz.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
3
Time step: 1 Sim. Time [au] 20.00
O -0.12083938E-02 0.11830991E+00 0.71734779E-03
H 0.76066994E+00 -0.47523035E+00 -0.76389825E-02
H -0.74248771E+00 -0.47367788E+00 -0.36354917E-02
3
Time step: 1 Sim. Time [au] 20.00
O 0.80985387E-04 0.11694942E+00 0.45556695E-03
H 0.76694201E+00 -0.47590015E+00 -0.74345295E-02
H -0.76900922E+00 -0.45157862E+00 0.31604714E-03
3
Time step: 1 Sim. Time [au] 20.00
O -0.33627926E-03 0.11839556E+00 0.65861546E-04
H 0.76654469E+00 -0.48039230E+00 -0.61359349E-02
H -0.76196819E+00 -0.46988778E+00 0.52038584E-02
3
Time step: 1 Sim. Time [au] 20.00
O 0.84435611E-03 0.11710864E+00 -0.66826162E-04
H 0.75418761E+00 -0.46371653E+00 -0.39391948E-02
H -0.76820777E+00 -0.46619053E+00 0.51147382E-02
3
Time step: 1 Sim. Time [au] 20.00
O 0.23959871E-02 0.11751167E+00 0.14205680E-02
H 0.75329894E+00 -0.46996488E+00 -0.73862035E-02
H -0.79179511E+00 -0.46635384E+00 -0.14960403E-01
3
Time step: 1 Sim. Time [au] 20.00
O -0.44172695E-03 0.11890421E+00 -0.11987099E-02
H 0.76595808E+00 -0.46702691E+00 0.10228920E-01
H -0.75989908E+00 -0.49126558E+00 0.88641280E-02
3
Time step: 1 Sim. Time [au] 20.00
O -0.30965923E-03 0.11754233E+00 -0.82965765E-03
H 0.76003487E+00 -0.46651093E+00 0.13038239E-01
H -0.75606282E+00 -0.47027618E+00 0.19951548E-03
3
Time step: 1 Sim. Time [au] 20.00
O 0.97105547E-03 0.11849638E+00 -0.14038561E-02
H 0.74873717E+00 -0.48604620E+00 0.12547242E-01
H -0.76507089E+00 -0.46581552E+00 0.97915471E-02
3
Time step: 1 Sim. Time [au] 20.00
O 0.20802456E-03 0.11737082E+00 -0.41768812E-03
H 0.75600119E+00 -0.46957273E+00 0.78963913E-02
H -0.76022122E+00 -0.46447187E+00 -0.11991736E-02
3
Time step: 1 Sim. Time [au] 20.00
O 0.98697769E-03 0.11894732E+00 -0.90415747E-03
H 0.74609216E+00 -0.47761489E+00 0.91261648E-02
H -0.76265778E+00 -0.48140610E+00 0.52847855E-02
Loading

0 comments on commit a8da962

Please sign in to comment.