Skip to content

Commit

Permalink
Functionalities to incorporate rotational averaging in i-PI (i-pi#372)
Browse files Browse the repository at this point in the history
* Added functions to manipulate rotations
* Added some utilities for rotational averaging
* GenericCell object handling non-triangular cells
This in no way allows using generic cells in i-PI, but allows us to handle arbitrary rotations from the forcefield to the socket
* Added an example, and options, to use lebedev grids in o3-averaging
* Implement defining and storing a separate PRNG for the random rotations
* smaller box size for geop examples
these make the tests run much faster
  • Loading branch information
ceriottm authored Dec 24, 2024
1 parent e0d5b7b commit 9f2c98d
Show file tree
Hide file tree
Showing 51 changed files with 979 additions and 48 deletions.
13 changes: 6 additions & 7 deletions drivers/f90/LJ.f90
Original file line number Diff line number Diff line change
Expand Up @@ -190,7 +190,7 @@ SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list,
pot = pot + pot_ij
DO k = 1, 3
DO l = k, 3
! Only the upper triangular elements calculated.
! only upper triangular part is returned
virial(k,l) = virial(k,l) + fij(k)*rij(l)
ENDDO
ENDDO
Expand All @@ -199,8 +199,8 @@ SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list,
start = index_list(i) + 1
ENDDO

! Assuming an upper-triangular vector matrix for the simulation box.
volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
! Works with a generic cell, even if usually it'll be upper-triuangular
volume = CELL_VOLUME(cell_h) !cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
CALL LJ_longrange(rc, sigma, eps, natoms, volume, pot_lr, vir_lr)
pot = pot + pot_lr
DO k = 1, 3
Expand Down Expand Up @@ -272,8 +272,7 @@ SUBROUTINE LJMix_getall(n_type2, rc, sigma, eps, natoms, atoms, cell_h, cell_ih,
forces(n_list(j),:) = forces(n_list(j),:) - fij
pot = pot + pot_ij
DO k = 1, 3
DO l = k, 3
! Only the upper triangular elements calculated.
DO l = k, 3 ! only upper-tri part is returned
virial(k,l) = virial(k,l) + fij(k)*rij(l)
ENDDO
ENDDO
Expand All @@ -282,8 +281,8 @@ SUBROUTINE LJMix_getall(n_type2, rc, sigma, eps, natoms, atoms, cell_h, cell_ih,
start = index_list(i) + 1
ENDDO

! Assuming an upper-triangular vector matrix for the simulation box.
volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
! Works with a generic cell, even if usually it'll be upper-triuangular
volume = CELL_VOLUME(cell_h) ! cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
CALL LJ_longrange(rc, sigma*(1-n_type2*0.4/natoms), eps*(1-n_type2*0.4/natoms), natoms, volume, pot_lr, vir_lr)
pot = pot + pot_lr
DO k = 1, 3
Expand Down
11 changes: 11 additions & 0 deletions drivers/f90/distance.f90
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,17 @@ MODULE DISTANCE

CONTAINS

DOUBLE PRECISION function CELL_VOLUME(m)

DOUBLE PRECISION m(3, 3)

CELL_VOLUME = DABS(m(1, 1) * (m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3)) - &
m(1, 2) * (m(2, 1) * m(3, 3) - m(3, 1) * m(2, 3)) + &
m(1, 3) * (m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2)) )

return
end

SUBROUTINE vector_separation(cell_h, cell_ih, ri, rj, rij, r2)
! Calculates the vector separating two atoms.
!
Expand Down
48 changes: 43 additions & 5 deletions drivers/f90/driver.f90
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ PROGRAM DRIVER
USE SG
USE PSWATER
USE F90SOCKETS, ONLY : open_socket, writebuffer, readbuffer, f_sleep
USE DISTANCE, only: CELL_VOLUME
IMPLICIT NONE

! SOCKET VARIABLES
Expand Down Expand Up @@ -191,13 +192,15 @@ PROGRAM DRIVER
vstyle = 64
ELSEIF (trim(cmdbuffer) == "qtip4pf-c-json-delta") THEN
vstyle = 65
ELSEIF (trim(cmdbuffer) == "noo3-h2o") THEN
vstyle = 70
ELSEIF (trim(cmdbuffer) == "gas") THEN
vstyle = 0 ! ideal gas
ELSEIF (trim(cmdbuffer) == "dummy") THEN
vstyle = 99 ! returns non-zero but otherwise meaningless values
ELSE
WRITE(*,*) " Unrecognized potential type ", trim(cmdbuffer)
WRITE(*,*) " Use -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4pf-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta] "
WRITE(*,*) " Use -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4pf-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta|noo3-h2o] "
STOP "ENDED"
ENDIF
ELSEIF (ccmd == 4) THEN
Expand Down Expand Up @@ -385,7 +388,7 @@ PROGRAM DRIVER
enddo
ENDIF
isinit = .true.
ELSEIF (vstyle == 1) THEN
ELSEIF (1 == vstyle) THEN
IF (par_count /= 3) THEN
WRITE(*,*) "Error: parameters not initialized correctly."
WRITE(*,*) "For LJ potential use -o sigma,epsilon,cutoff "
Expand Down Expand Up @@ -519,8 +522,8 @@ PROGRAM DRIVER
! units and storage mode used in the driver.
cell_h = transpose(cell_h)
cell_ih = transpose(cell_ih)
! We assume an upper triangular cell-vector matrix
volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
! We compute for a generic cell, just in case (even though usually i-PI passes an upper triangular cell-vector matrix)
volume = CELL_VOLUME(cell_h) !cell_h(1,1)*cell_h(2,2)*cell_h(3,3)

CALL readbuffer(socket, cbuf) ! The number of atoms in the cell
IF (verbose > 1) WRITE(*,*) " !read!=> cbuf: ", cbuf
Expand Down Expand Up @@ -790,6 +793,41 @@ PROGRAM DRIVER
CALL LJ_getall(rc, 2.1d0, -12d-6, nat, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
CALL LJ_getall(rc, 2.5d0, 1d-6, nat, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
ENDIF
ELSEIF (70 == vstyle) THEN
! potential that can be applied on top of water molecules to add a orientation-dependent term to
! test non-equivariant terms in the potential
vpars(1) = cell_h(1,1)
vpars(2) = cell_h(2,2)
vpars(3) = cell_h(3,3)

pot = 0
forces = 0
DO i=1,nat,3
dip = 0.5 * (atoms(i+1,:) + atoms(i+2,:)) - atoms(i,:)
vpars(4) = sqrt(dip(1)**2+dip(2)**2+dip(3)**2)
vpars(1) = exp(dip(1)/vpars(4)*1.0)
vpars(2) = exp(dip(2)/vpars(4)*0.5)
vpars(3) = exp(dip(3)/vpars(4)*2.0)
pot = pot + 2e-2*(vpars(1)+vpars(2)+vpars(3))
vpars(1) = 2e-2*vpars(1)*1.0/vpars(4)**3
vpars(2) = 2e-2*vpars(2)*0.5/vpars(4)**3
vpars(3) = 2e-2*vpars(3)*2.0/vpars(4)**3
! gradients on O
vpars(4) = -vpars(1)*(dip(2)**2+dip(3)**2) + vpars(2)*dip(1)*dip(2) + vpars(3)*dip(3)*dip(1)
vpars(5) = -vpars(2)*(dip(1)**2+dip(3)**2) + vpars(1)*dip(1)*dip(2) + vpars(3)*dip(2)*dip(3)
vpars(6) = -vpars(3)*(dip(1)**2+dip(2)**2) + vpars(1)*dip(1)*dip(3) + vpars(2)*dip(3)*dip(2)
forces(i,1) = forces(i,1) - vpars(4)
forces(i,2) = forces(i,2) - vpars(5)
forces(i,3) = forces(i,3) - vpars(6)
! gradients from H are just from Newton's law
forces(i+1,1) = forces(i+1,1) + vpars(4)*0.5
forces(i+1,2) = forces(i+1,2) + vpars(5)*0.5
forces(i+1,3) = forces(i+1,3) + vpars(6)*0.5
forces(i+2,1) = forces(i+2,1) + vpars(4)*0.5
forces(i+2,2) = forces(i+2,2) + vpars(5)*0.5
forces(i+2,3) = forces(i+2,3) + vpars(6)*0.5
ENDDO

ELSEIF (vstyle == 11) THEN ! efield potential.
IF (mod(nat,3)/=0) THEN
WRITE(*,*) " Expecting water molecules O H H O H H O H H but got ", nat, "atoms"
Expand Down Expand Up @@ -1055,7 +1093,7 @@ PROGRAM DRIVER
SUBROUTINE helpmessage
! Help banner

WRITE(*,*) " SYNTAX: driver.x [-u] -a address [-p port] -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4p-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta]"
WRITE(*,*) " SYNTAX: driver.x [-u] -a address [-p port] -m [dummy|gas|lj|sg|harm|harm3d|morse|morsedia|zundel|qtip4pf|pswater|lepsm1|lepsm2|qtip4p-efield|eckart|ch4hcbe|ljpolymer|MB|doublewell|doublewell_1D|water_dip_pol|harmonic_bath|meanfield_bath|ljmix|qtip4pf-sr|qtip4pf-c-1|qtip4pf-c-2|qtip4pf-c-json|qtip4pf-c-1-delta|qtip4pf-c-2-delta|qtip4pf-c-json-delta|noo3-h2o]"
WRITE(*,*) " -o 'comma_separated_parameters' [-S sockets_prefix] [-v] "
WRITE(*,*) ""
WRITE(*,*) " For LJ potential use -o sigma,epsilon,cutoff "
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
<prng>
<seed>31415</seed>
</prng>
<ffcommittee name="committee">
<ffcommittee name="committee" threaded="False">
<alpha> 1.0 </alpha>
<active_output> active </active_output>
<active_thresh units="electronvolt"> 0.05 </active_thresh>
Expand Down
21 changes: 21 additions & 0 deletions examples/features/o3_averaging/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
O3 rotational averaging.
========================

Examples of the averaging module to handle non-equivariant potentials.
The `<ffrotations>` socket allows to evaluate a potential multiple times,
using random rotations or a grid of rotations, to reduce the impact of the
lack of rotational equivariance on the quality of the dynamics.
The examples use the `noo3_h2o` potential, that assumes the structure is
a sequence of H2O molecules and adds a orientation-dependent potential to each
water molecule.


`noo3`: no correction is applied. The structure of the liquid will be badly
disrupted.
`random`: the structure is randomly oriented at each evaluation.
This effectively introduces a noise in the forces, that must be contrasted
by an aggressive thermostat. There is a large drift in the conserved quantity
but the mean potential is OK and water molecules do not orient artificially
`grid`: the potential is averaged over a grid of rotations.
This reduces the orientational error, but comes at the cost of many additional
energy evaluations per step.
42 changes: 42 additions & 0 deletions examples/features/o3_averaging/direct_liquid_water/input.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
<simulation verbosity='medium'>
<output prefix='simulation'>
<properties stride='1' filename='out'> [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] </properties>
<trajectory filename='noo3_pots' stride='2' cell_units='angstrom' extra_type="o3grid_pots"> extras_component_raw(1) </trajectory>
<trajectory filename='pos' stride='20' cell_units='angstrom'> positions{angstrom} </trajectory>
<checkpoint stride='200'/>
</output>
<total_steps>1000</total_steps>
<prng>
<seed>31415</seed>
</prng>
<ffsocket name='driver-qtip4pf' mode='unix' pbc='false'>
<address>h2o-base</address>
</ffsocket>
<ffrotations name='driver-noo3' pbc='false'>
<ffdirect name='base-noo3'>
<pes> dummy </pes> <!-- This is just a dummy ff -->
</ffdirect>
<random> True </random>
</ffrotations>
<system>
<initialize nbeads='1'>
<file mode='xyz'> water_216.xyz </file>
<velocities mode='thermal' units='kelvin'> 300 </velocities>
</initialize>
<forces>
<force forcefield='driver-qtip4pf'/>
<force forcefield='driver-noo3' weight="1"/>
</forces>
<motion mode='dynamics'>
<dynamics mode='nvt'>
<timestep units='femtosecond'> 0.5 </timestep>
<thermostat mode='langevin'>
<tau units='femtosecond'> 10 </tau>
</thermostat>
</dynamics>
</motion>
<ensemble>
<temperature units='kelvin'> 300 </temperature>
</ensemble>
</system>
</simulation>
16 changes: 16 additions & 0 deletions examples/features/o3_averaging/direct_liquid_water/run.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
ipi=i-pi
driver1="i-pi-driver -m qtip4pf -u -a h2o-base"
sleep_time=4

${ipi} input.xml > log.i-pi &
echo "# i-PI is running"

echo "# Waiting for ${sleep_time} (s) before executing driver"
sleep ${sleep_time}

${driver1} > /dev/null &
echo "# Driver is running"

wait

echo "# Simulation complete"
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
#Settings file for automatic test
driver_model qtip4pf
nsteps 10
44 changes: 44 additions & 0 deletions examples/features/o3_averaging/grid_liquid_water/input.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
<simulation verbosity='medium'>
<output prefix='simulation'>
<properties stride='1' filename='out'> [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] </properties>
<trajectory filename='noo3_pots' stride='2' cell_units='angstrom' extra_type="o3grid_pots"> extras_component_raw(1) </trajectory>
<trajectory filename='pos' stride='20' cell_units='angstrom'> positions{angstrom} </trajectory>
<checkpoint stride='200'/>
</output>
<total_steps>1000</total_steps>
<prng>
<seed>31415</seed>
</prng>
<ffsocket name='driver-qtip4pf' mode='unix' pbc='false'>
<address>h2o-base</address>
</ffsocket>
<ffrotations name='driver-noo3' pbc='false'>
<ffsocket name='base-noo3' mode='unix' pbc='false'>
<address>h2o-noo3</address>
</ffsocket>
<grid_order> 2 </grid_order>
<grid_mode> legendre </grid_mode>
<inversion> True </inversion>
</ffrotations>
<system>
<initialize nbeads='1'>
<file mode='xyz'> water_216.xyz </file>
<velocities mode='thermal' units='kelvin'> 300 </velocities>
</initialize>
<forces>
<force forcefield='driver-qtip4pf'/>
<force forcefield='driver-noo3' weight="1.0"/>
</forces>
<motion mode='dynamics'>
<dynamics mode='nvt'>
<timestep units='femtosecond'> 0.5 </timestep>
<thermostat mode='langevin'>
<tau units='femtosecond'> 100 </tau>
</thermostat>
</dynamics>
</motion>
<ensemble>
<temperature units='kelvin'> 300 </temperature>
</ensemble>
</system>
</simulation>
18 changes: 18 additions & 0 deletions examples/features/o3_averaging/grid_liquid_water/run.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
ipi=i-pi
driver1="i-pi-driver -m noo3-h2o -u -a h2o-noo3"
driver2="i-pi-driver -m qtip4pf -u -a h2o-base"
sleep_time=4

${ipi} input.xml > log.i-pi &
echo "# i-PI is running"

echo "# Waiting for ${sleep_time} (s) before executing driver"
sleep ${sleep_time}

${driver1} > /dev/null &
${driver2} > /dev/null &
echo "# Driver is running"

wait

echo "# Simulation complete"
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
#Settings file for automatic test
driver_model qtip4pf
driver_model2 noo3-h2o
nsteps 10
43 changes: 43 additions & 0 deletions examples/features/o3_averaging/lebedev_liquid_water/input.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
<simulation verbosity='high'>
<output prefix='simulation'>
<properties stride='1' filename='out'> [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] </properties>
<trajectory filename='noo3_pots' stride='2' cell_units='angstrom' extra_type="o3grid_pots"> extras_component_raw(1) </trajectory>
<trajectory filename='pos' stride='20' cell_units='angstrom'> positions{angstrom} </trajectory>
<checkpoint stride='200'/>
</output>
<total_steps>1000</total_steps>
<prng>
<seed>31415</seed>
</prng>
<ffsocket name='driver-qtip4pf' mode='unix' pbc='false'>
<address>h2o-base</address>
</ffsocket>
<ffrotations name='driver-noo3' pbc='false'>
<ffsocket name='base-noo3' mode='unix' pbc='false'>
<address>h2o-noo3</address>
</ffsocket>
<grid_order> 3 </grid_order> <grid_mode> lebedev </grid_mode>
<inversion> True </inversion>
</ffrotations>
<system>
<initialize nbeads='1'>
<file mode='xyz'> water_216.xyz </file>
<velocities mode='thermal' units='kelvin'> 300 </velocities>
</initialize>
<forces>
<force forcefield='driver-qtip4pf'/>
<force forcefield='driver-noo3' weight="1.0"/>
</forces>
<motion mode='dynamics'>
<dynamics mode='nvt'>
<timestep units='femtosecond'> 0.5 </timestep>
<thermostat mode='langevin'>
<tau units='femtosecond'> 100 </tau>
</thermostat>
</dynamics>
</motion>
<ensemble>
<temperature units='kelvin'> 300 </temperature>
</ensemble>
</system>
</simulation>
18 changes: 18 additions & 0 deletions examples/features/o3_averaging/lebedev_liquid_water/run.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
ipi=i-pi
driver1="i-pi-driver -m noo3-h2o -u -a h2o-noo3"
driver2="i-pi-driver -m qtip4pf -u -a h2o-base"
sleep_time=4

${ipi} input.xml > log.i-pi &
echo "# i-PI is running"

echo "# Waiting for ${sleep_time} (s) before executing driver"
sleep ${sleep_time}

${driver1} > /dev/null &
${driver2} > /dev/null &
echo "# Driver is running"

wait

echo "# Simulation complete"
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
#Settings file for automatic test
driver_model qtip4pf
driver_model2 noo3-h2o
nsteps 10
Loading

0 comments on commit 9f2c98d

Please sign in to comment.