Skip to content

Commit

Permalink
Implement numerical forces for H2O Schwenke potential (#174)
Browse files Browse the repository at this point in the history 
* Add h2o_scan OH bond scans for various potentials

---------

Co-authored-by: Daniel Hollas <daniel.hollas@bristol.ac.uk>
  • Loading branch information
@mbarnfield63 @danielhollas
mbarnfield63 and danielhollas authored Mar 11, 2024
1 parent 2b6a0b3 commit b5383db
Show file tree
Hide file tree
Showing 18 changed files with 1,870 additions and 8 deletions.
68 changes: 63 additions & 5 deletions src/force_h2o.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,6 @@ subroutine force_h2o_schwenke(x, y, z, fx, fy, fz, Eclas, natom, nbeads)

! For Path Integrals, the final energy of the PI necklace
! is an average over all beads.
! TODO: Forces need to be appropriately scaled as well!
do iw = 1, nbeads
Eclas = Eclas + Epot(iw)
end do
Expand Down Expand Up @@ -135,19 +134,78 @@ end subroutine force_h2o_cvrqd

! TODO: Implement numerical forces generally for all potentials
! For now, they can be implemented here and hardcoded for a specific H2O potential
subroutine numerical_forces(x, y, z, Epot, fx, fy, fz, natom, nbeads)
subroutine numerical_forces(x, y, z, fx, fy, fz, Epot, natom, nbeads)
real(DP), intent(in) :: x(natom, nbeads)
real(DP), intent(in) :: y(natom, nbeads)
real(DP), intent(in) :: z(natom, nbeads)
real(DP), intent(inout) :: fx(natom, nbeads)
real(DP), intent(inout) :: fy(natom, nbeads)
real(DP), intent(inout) :: fz(natom, nbeads)
integer, intent(in) :: natom, nbeads

! Create new copies of arrays
real(DP) :: x_new_forward(natom, nbeads)
real(DP) :: y_new_forward(natom, nbeads)
real(DP) :: z_new_forward(natom, nbeads)

! This is the energy for the currrent geometry that has already been calculated
real(DP), intent(in) :: Epot(nbeads)
integer, intent(in) :: natom, nbeads

! Need to implement numerical forces first
call fatal_error(__FILE__, __LINE__, 'Numerical forces not yet implemented!')
! Internal water coordinates
real(DP) :: new_rOH1, new_rOH2, new_aHOH_rad
real(DP) :: new_rij(1, 3)

! Schwenke calculated peterbed geometry energy
real(DP) :: Epot_delta(1)

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

! Calculate forces numerically using central differences
do j = 1, nbeads

! Save the original energy
Eclas_orig = Epot(j)

do i = 1, natom

do k = 1, 3 ! x, y, z

! Copy the original atom coordinates
x_new_forward(:, :) = x(:, :)
y_new_forward(:, :) = y(:, :)
z_new_forward(:, :) = z(:, :)

! Move the atom forwards
select case (k)
case (1)
x_new_forward(i, j) = x_new_forward(i, j) + delta
case (2)
y_new_forward(i, j) = y_new_forward(i, j) + delta
case (3)
z_new_forward(i, j) = z_new_forward(i, j) + delta
end select

! Calculate the energy for the forward perturbed geometry
call get_internal_coords(x_new_forward, y_new_forward, z_new_forward, j, new_rOH1, new_rOH2, new_aHOH_rad)

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

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

! Calculate the numerical force
select case (k)
case (1)
fx(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
case (2)
fy(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
case (3)
fz(i, j) = -(Epot_delta(1) - Eclas_orig) / delta
end select
end do
end do
end do
end subroutine numerical_forces

end module mod_force_h2o
67 changes: 67 additions & 0 deletions tests/H2O_CVRQD/h2o_scan_cvrqd.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
! To compile:
! gfortran -c h2o_scan_cvrqd.f90
! gfortran -fopenmp -L./ -L../water_potentials/ libabin.a ../water_potentials/libwater.a h2o_scan_cvrqd.o -labin -lwater -lm -lstdc++ -o cvrqd_scan
program cvrqd_scan
use mod_force_h2o
use mod_const, only: DP, ANG
use mod_init, only: read_xyz_file
implicit none

integer, parameter :: NATOM = 3
integer, parameter :: NBEADS = 1
integer, parameter :: WATPOT = 1

! Displacement for H atom at each timestep
real(DP) :: delta = 0.0005_DP
integer, parameter :: NSTEPS = 5000

real(DP) :: x(3, 1), y(3, 1), z(3, 1)
real(DP) :: fx(3, 1), fy(3, 1), fz(3, 1)
real(DP) :: energy, rOH
integer :: i, funit
character(len=2) :: atom_names(3)
character(len=256) :: fname
!character(len=20) :: h2opot = 'cvrqd'

atom_names = (/"O", "H", "H"/)

fname = "water.xyz"
open (newunit=funit, file=fname, action="read")
call read_xyz_file(funit, fname, atom_names, natom, nbeads, x, y, z)
close (funit)

! Convert from Angstromgs to Bohr (atomic units)
x(:, 1) = x(:, 1) * ANG
y(:, 1) = y(:, 1) * ANG
z(:, 1) = z(:, 1) * ANG

! print for output
print*,"Step ", "rOH ", "Energy"

open (99, file="water_geometries.xyz", action="write")
! steps for stretching OH bond
do i = 1, NSTEPS
! Calculate the initial OH bond length (rOH) from initial bond positions
rOH = sqrt((x(2, 1) - x(1, 1))**2 + (y(2, 1) - y(1, 1))**2 + (z(2, 1) - z(1, 1))**2)

! Modify the OH bond length by changing the position of a single hydrogen atom
x(2, 1) = x(2, 1) + delta

! Call the force_water function to calculate forces and energy
energy = 0.0D0
call force_h2o(x, y, z, fx, fy, fz, energy, natom, nbeads)

! Print the current OH bond length and energy
print '(I4,F16.8,E18.8)', i, rOH, energy

! Print XYZ geometry
write (99, '(I1)') NATOM
write (99, *)
write (99, '(A,F12.6,F12.6,F12.6)') "O ", x(1, 1) / ANG, y(1, 1) / ANG, z(1, 1) / ANG
write (99, '(A,F12.6,F12.6,F12.6)') "H ", x(2, 1) / ANG, y(2, 1) / ANG, z(2, 1) / ANG
write (99, '(A,F12.6,F12.6,F12.6)') "H ", x(3, 1) / ANG, y(3, 1) / ANG, z(3, 1) / ANG
end do

! close file
close (99)
end program
1 change: 0 additions & 1 deletion tests/H2O_SCHWENKE/ERROR.ref
View file

This file was deleted.

2 changes: 2 additions & 0 deletions tests/H2O_SCHWENKE/energies.dat.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# Time[fs] E-potential E-kinetic E-Total E-Total-Avg
0.48 0.6381703173E-05 0.7149142696E-06 0.7096617443E-05 0.7096617443E-05
5 changes: 5 additions & 0 deletions tests/H2O_SCHWENKE/forces.xyz.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
3
net force: -0.35413E-04 -0.22257E-04 -0.10937E-06 torque force: -0.29376E-07 0.13497E-10 0.12718E-04
O -0.1752322176E-04 -0.2737840439E-02 -0.6114886375E-07
H -0.1019196611E-02 0.1357557850E-02 -0.2410618973E-07
H 0.1001307070E-02 0.1358025502E-02 -0.2411562662E-07
67 changes: 67 additions & 0 deletions tests/H2O_SCHWENKE/h2o_scan_qtip4p.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
! To compile:
! gfortran -c h2o_scan.f90
! gfortran -fopenmp -L./ -L../water_potentials/ libabin.a ../water_potentials/libwater.a h2o_scan.o -labin -lwater -lm -lstdc++ -o scan
program scan
!use mod_water
use mod_const, only: DP, ANG
use mod_init, only: read_xyz_file
implicit none
external :: force_water

integer, parameter :: NATOM = 3
integer, parameter :: NBEADS = 1
integer, parameter :: WATPOT = 1

! Displacement for H atom at each timestep
real(DP) :: delta = 0.0005_DP
integer, parameter :: NSTEPS = 5000

real(DP) :: x(3, 1), y(3, 1), z(3, 1)
real(DP) :: fx(3, 1), fy(3, 1), fz(3, 1)
real(DP) :: energy, rOH
integer :: i, funit
character(len=2) :: atom_names(3)
character(len=256) :: fname

atom_names = (/"O", "H", "H"/)

fname = "water.xyz"
open (newunit=funit, file=fname, action="read")
call read_xyz_file(funit, fname, atom_names, natom, nbeads, x, y, z)
close (funit)

! Convert from Angstromgs to Bohr (atomic units)
x(:, 1) = x(:, 1) * ANG
y(:, 1) = y(:, 1) * ANG
z(:, 1) = z(:, 1) * ANG

! print for output
print*,"Step ", "rOH ", "Energy"

open (99, file="water_geometries.xyz", action="write")
! steps for stretching OH bond
do i = 1, NSTEPS
! Calculate the initial OH bond length (rOH) from initial bond positions
rOH = sqrt((x(2, 1) - x(1, 1))**2 + (y(2, 1) - y(1, 1))**2 + (z(2, 1) - z(1, 1))**2)

! Modify the OH bond length by changing the position of a single hydrogen atom
x(2, 1) = x(2, 1) + delta

! Call the force_water function to calculate forces and energy
energy = 0.0D0
call force_water(x, y, z, fx, fy, fz, energy, natom, nbeads, watpot)

! Print the current OH bond length and energy
print '(I4,F16.8,E18.8)', i, rOH, energy

! Print XYZ geometry
write (99, '(I1)') NATOM
write (99, *)
write (99, '(A,F12.6,F12.6,F12.6)') "O ", x(1, 1) / ANG, y(1, 1) / ANG, z(1, 1) / ANG
write (99, '(A,F12.6,F12.6,F12.6)') "H ", x(2, 1) / ANG, y(2, 1) / ANG, z(2, 1) / ANG
write (99, '(A,F12.6,F12.6,F12.6)') "H ", x(3, 1) / ANG, y(3, 1) / ANG, z(3, 1) / ANG
end do

! close file
close (99)
end program
66 changes: 66 additions & 0 deletions tests/H2O_SCHWENKE/h2o_scan_schwenke.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
! To compile:
! gfortran -c h2o_scan_schwenke.f90
! gfortran -fopenmp -L./ -L../water_potentials/ libabin.a ../water_potentials/libwater.a h2o_scan_schwenke.o -labin -lwater -lm -lstdc++ -o schwenke_scan
program schwenke_scan
use mod_force_h2o
use mod_const, only: DP, ANG
use mod_init, only: read_xyz_file
implicit none

integer, parameter :: NATOM = 3
integer, parameter :: NBEADS = 1
integer, parameter :: WATPOT = 1

! Displacement for H atom at each timestep
real(DP) :: delta = 0.0005_DP
integer, parameter :: NSTEPS = 5000

real(DP) :: x(3, 1), y(3, 1), z(3, 1)
real(DP) :: fx(3, 1), fy(3, 1), fz(3, 1)
real(DP) :: energy, rOH
integer :: i, funit
character(len=2) :: atom_names(3)
character(len=256) :: fname

atom_names = (/"O", "H", "H"/)

fname = "water.xyz"
open (newunit=funit, file=fname, action="read")
call read_xyz_file(funit, fname, atom_names, natom, nbeads, x, y, z)
close (funit)

! Convert from Angstromgs to Bohr (atomic units)
x(:, 1) = x(:, 1) * ANG
y(:, 1) = y(:, 1) * ANG
z(:, 1) = z(:, 1) * ANG

! print for output
print*,"Step ", "rOH ", "Energy"

open (99, file="water_geometries.xyz", action="write")
! steps for stretching OH bond
do i = 1, NSTEPS
! Calculate the initial OH bond length (rOH) from initial bond positions
rOH = sqrt((x(2, 1) - x(1, 1))**2 + (y(2, 1) - y(1, 1))**2 + (z(2, 1) - z(1, 1))**2)

! Modify the OH bond length by changing the position of a single hydrogen atom
x(2, 1) = x(2, 1) + delta

! Call the force_water function to calculate forces and energy
energy = 0.0D0
call force_h2o(x, y, z, fx, fy, fz, energy, natom, nbeads)

! Print the current OH bond length and energy
print '(I4,F16.8,E18.8)', i, rOH, energy

! Print XYZ geometry
write (99, '(I1)') NATOM
write (99, *)
write (99, '(A,F12.6,F12.6,F12.6)') "O ", x(1, 1) / ANG, y(1, 1) / ANG, z(1, 1) / ANG
write (99, '(A,F12.6,F12.6,F12.6)') "H ", x(2, 1) / ANG, y(2, 1) / ANG, z(2, 1) / ANG
write (99, '(A,F12.6,F12.6,F12.6)') "H ", x(3, 1) / ANG, y(3, 1) / ANG, z(3, 1) / ANG
end do

! close file
close (99)
end program
3 changes: 2 additions & 1 deletion tests/H2O_SCHWENKE/input.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,11 @@ h2opot='schwenke'
mdtype='MD', ! classical MD
dt=20., ! number of steps and timestep
nstep=1
nwritef=1, ! write forces
/

&nhcopt
inose=0, ! Thermostating: Nose-Hoover 1, microcanonical 0,GLE 2, LE 3
temp=100
temp=0
rem_comrot=.true. ! this is a default value, remove rotations at the beginning
/
5 changes: 5 additions & 0 deletions tests/H2O_SCHWENKE/movie.xyz.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
3
Time step: 1 Sim. Time [au] 20.00
O -0.61360661E-07 0.11798960E+00 -0.23219057E-09
H 0.75693652E+00 -0.47191813E+00 -0.14713543E-08
H -0.75693758E+00 -0.47191813E+00 -0.14713543E-08
2 changes: 2 additions & 0 deletions tests/H2O_SCHWENKE/temper.dat.ref
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# Time[fs] Temperature T-Average Conserved_quantity_of_thermostat
0.48 0.05 0.05
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
Loading

0 comments on commit b5383db

Please sign in to comment.