diff --git a/drivers/f90/LJ.f90 b/drivers/f90/LJ.f90 index 67343ef08..ab21d7b13 100644 --- a/drivers/f90/LJ.f90 +++ b/drivers/f90/LJ.f90 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/drivers/f90/distance.f90 b/drivers/f90/distance.f90 index c15213030..c74c2a624 100644 --- a/drivers/f90/distance.f90 +++ b/drivers/f90/distance.f90 @@ -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. ! diff --git a/drivers/f90/driver.f90 b/drivers/f90/driver.f90 index c1554f0fd..5f210c6da 100644 --- a/drivers/f90/driver.f90 +++ b/drivers/f90/driver.f90 @@ -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 @@ -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 @@ -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 " @@ -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 @@ -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" @@ -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 " diff --git a/examples/features/committee_models/thermodynamic_uq/input.xml b/examples/features/committee_models/thermodynamic_uq/input.xml index ef96ac06d..df838e593 100644 --- a/examples/features/committee_models/thermodynamic_uq/input.xml +++ b/examples/features/committee_models/thermodynamic_uq/input.xml @@ -12,7 +12,7 @@ 31415 - + 1.0 active 0.05 diff --git a/examples/features/o3_averaging/README.md b/examples/features/o3_averaging/README.md new file mode 100644 index 000000000..85cbd576c --- /dev/null +++ b/examples/features/o3_averaging/README.md @@ -0,0 +1,21 @@ +O3 rotational averaging. +======================== + +Examples of the averaging module to handle non-equivariant potentials. +The `` 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. diff --git a/examples/features/o3_averaging/direct_liquid_water/input.xml b/examples/features/o3_averaging/direct_liquid_water/input.xml new file mode 100644 index 000000000..2d85ee633 --- /dev/null +++ b/examples/features/o3_averaging/direct_liquid_water/input.xml @@ -0,0 +1,42 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] + extras_component_raw(1) + positions{angstrom} + + + 1000 + + 31415 + + +
h2o-base
+ + + + dummy + + True + + + + water_216.xyz + 300 + + + + + + + + 0.5 + + 10 + + + + + 300 + + + diff --git a/examples/features/o3_averaging/direct_liquid_water/run.sh b/examples/features/o3_averaging/direct_liquid_water/run.sh new file mode 100644 index 000000000..be4dcec32 --- /dev/null +++ b/examples/features/o3_averaging/direct_liquid_water/run.sh @@ -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" diff --git a/examples/features/o3_averaging/direct_liquid_water/test_settings.dat b/examples/features/o3_averaging/direct_liquid_water/test_settings.dat new file mode 100644 index 000000000..931f44954 --- /dev/null +++ b/examples/features/o3_averaging/direct_liquid_water/test_settings.dat @@ -0,0 +1,3 @@ +#Settings file for automatic test +driver_model qtip4pf +nsteps 10 diff --git a/examples/features/o3_averaging/direct_liquid_water/water_216.xyz b/examples/features/o3_averaging/direct_liquid_water/water_216.xyz new file mode 120000 index 000000000..58a9ede12 --- /dev/null +++ b/examples/features/o3_averaging/direct_liquid_water/water_216.xyz @@ -0,0 +1 @@ +../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/features/o3_averaging/grid_liquid_water/input.xml b/examples/features/o3_averaging/grid_liquid_water/input.xml new file mode 100644 index 000000000..fe86ed9e5 --- /dev/null +++ b/examples/features/o3_averaging/grid_liquid_water/input.xml @@ -0,0 +1,44 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] + extras_component_raw(1) + positions{angstrom} + + + 1000 + + 31415 + + +
h2o-base
+ + + +
h2o-noo3
+ + 2 + legendre + True + + + + water_216.xyz + 300 + + + + + + + + 0.5 + + 100 + + + + + 300 + + + diff --git a/examples/features/o3_averaging/grid_liquid_water/run.sh b/examples/features/o3_averaging/grid_liquid_water/run.sh new file mode 100644 index 000000000..504a1eb18 --- /dev/null +++ b/examples/features/o3_averaging/grid_liquid_water/run.sh @@ -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" diff --git a/examples/features/o3_averaging/grid_liquid_water/test_settings.dat b/examples/features/o3_averaging/grid_liquid_water/test_settings.dat new file mode 100644 index 000000000..5eb138751 --- /dev/null +++ b/examples/features/o3_averaging/grid_liquid_water/test_settings.dat @@ -0,0 +1,4 @@ +#Settings file for automatic test +driver_model qtip4pf +driver_model2 noo3-h2o +nsteps 10 diff --git a/examples/features/o3_averaging/grid_liquid_water/water_216.xyz b/examples/features/o3_averaging/grid_liquid_water/water_216.xyz new file mode 120000 index 000000000..58a9ede12 --- /dev/null +++ b/examples/features/o3_averaging/grid_liquid_water/water_216.xyz @@ -0,0 +1 @@ +../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/features/o3_averaging/lebedev_liquid_water/input.xml b/examples/features/o3_averaging/lebedev_liquid_water/input.xml new file mode 100644 index 000000000..118ffd54b --- /dev/null +++ b/examples/features/o3_averaging/lebedev_liquid_water/input.xml @@ -0,0 +1,43 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] + extras_component_raw(1) + positions{angstrom} + + + 1000 + + 31415 + + +
h2o-base
+ + + +
h2o-noo3
+ + 3 lebedev + True + + + + water_216.xyz + 300 + + + + + + + + 0.5 + + 100 + + + + + 300 + + + diff --git a/examples/features/o3_averaging/lebedev_liquid_water/run.sh b/examples/features/o3_averaging/lebedev_liquid_water/run.sh new file mode 100644 index 000000000..504a1eb18 --- /dev/null +++ b/examples/features/o3_averaging/lebedev_liquid_water/run.sh @@ -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" diff --git a/examples/features/o3_averaging/lebedev_liquid_water/test_settings.dat b/examples/features/o3_averaging/lebedev_liquid_water/test_settings.dat new file mode 100644 index 000000000..5eb138751 --- /dev/null +++ b/examples/features/o3_averaging/lebedev_liquid_water/test_settings.dat @@ -0,0 +1,4 @@ +#Settings file for automatic test +driver_model qtip4pf +driver_model2 noo3-h2o +nsteps 10 diff --git a/examples/features/o3_averaging/lebedev_liquid_water/water_216.xyz b/examples/features/o3_averaging/lebedev_liquid_water/water_216.xyz new file mode 120000 index 000000000..58a9ede12 --- /dev/null +++ b/examples/features/o3_averaging/lebedev_liquid_water/water_216.xyz @@ -0,0 +1 @@ +../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/features/o3_averaging/noo3_liquid_water/input.xml b/examples/features/o3_averaging/noo3_liquid_water/input.xml new file mode 100644 index 000000000..f6b509717 --- /dev/null +++ b/examples/features/o3_averaging/noo3_liquid_water/input.xml @@ -0,0 +1,38 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] + positions{angstrom} + + + 1000 + + 31415 + + +
h2o-base
+ + +
h2o-noo3
+ + + + water_216.xyz + 300 + + + + + + + + 0.5 + + 100 + + + + + 300 + + + diff --git a/examples/features/o3_averaging/noo3_liquid_water/run.sh b/examples/features/o3_averaging/noo3_liquid_water/run.sh new file mode 100644 index 000000000..504a1eb18 --- /dev/null +++ b/examples/features/o3_averaging/noo3_liquid_water/run.sh @@ -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" diff --git a/examples/features/o3_averaging/noo3_liquid_water/test_settings.dat b/examples/features/o3_averaging/noo3_liquid_water/test_settings.dat new file mode 100644 index 000000000..5eb138751 --- /dev/null +++ b/examples/features/o3_averaging/noo3_liquid_water/test_settings.dat @@ -0,0 +1,4 @@ +#Settings file for automatic test +driver_model qtip4pf +driver_model2 noo3-h2o +nsteps 10 diff --git a/examples/features/o3_averaging/noo3_liquid_water/water_216.xyz b/examples/features/o3_averaging/noo3_liquid_water/water_216.xyz new file mode 120000 index 000000000..58a9ede12 --- /dev/null +++ b/examples/features/o3_averaging/noo3_liquid_water/water_216.xyz @@ -0,0 +1 @@ +../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/features/o3_averaging/random_liquid_water/input.xml b/examples/features/o3_averaging/random_liquid_water/input.xml new file mode 100644 index 000000000..e8b196b7e --- /dev/null +++ b/examples/features/o3_averaging/random_liquid_water/input.xml @@ -0,0 +1,43 @@ + + + [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_md{electronvolt}, potential{electronvolt}, pressure_md{megapascal}, stress_md{megapascal} ] + extras_component_raw(1) + positions{angstrom} + + + 1000 + + 31415 + + +
h2o-base
+ + + +
h2o-noo3
+ + True + 54321 + + + + water_216.xyz + 300 + + + + + + + + 0.5 + + 10 + + + + + 300 + + + diff --git a/examples/features/o3_averaging/random_liquid_water/run.sh b/examples/features/o3_averaging/random_liquid_water/run.sh new file mode 100644 index 000000000..504a1eb18 --- /dev/null +++ b/examples/features/o3_averaging/random_liquid_water/run.sh @@ -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" diff --git a/examples/features/o3_averaging/random_liquid_water/test_settings.dat b/examples/features/o3_averaging/random_liquid_water/test_settings.dat new file mode 100644 index 000000000..5eb138751 --- /dev/null +++ b/examples/features/o3_averaging/random_liquid_water/test_settings.dat @@ -0,0 +1,4 @@ +#Settings file for automatic test +driver_model qtip4pf +driver_model2 noo3-h2o +nsteps 10 diff --git a/examples/features/o3_averaging/random_liquid_water/water_216.xyz b/examples/features/o3_averaging/random_liquid_water/water_216.xyz new file mode 120000 index 000000000..58a9ede12 --- /dev/null +++ b/examples/features/o3_averaging/random_liquid_water/water_216.xyz @@ -0,0 +1 @@ +../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/temp/lammps/h2o-piglet.2/data.lmp b/examples/temp/lammps/h2o-piglet.2/data.lmp index 13c75e993..f91ddfb1d 100644 --- a/examples/temp/lammps/h2o-piglet.2/data.lmp +++ b/examples/temp/lammps/h2o-piglet.2/data.lmp @@ -11,6 +11,7 @@ LAMMPS Description 0 35.233 xlo xhi 0 35.233 ylo yhi 0 35.233 zlo zhi + 1e-4 1e-4 1e-4 xy xz yz Masses diff --git a/examples/temp/pes/geop/waterbox/bfgs/input.xml b/examples/temp/pes/geop/waterbox/bfgs/input.xml index 9e9c3897d..643d27e13 100644 --- a/examples/temp/pes/geop/waterbox/bfgs/input.xml +++ b/examples/temp/pes/geop/waterbox/bfgs/input.xml @@ -12,7 +12,7 @@ - water_216.xyz + water_64.xyz diff --git a/examples/temp/pes/geop/waterbox/bfgs/water_216.xyz b/examples/temp/pes/geop/waterbox/bfgs/water_216.xyz deleted file mode 120000 index 95291df24..000000000 --- a/examples/temp/pes/geop/waterbox/bfgs/water_216.xyz +++ /dev/null @@ -1 +0,0 @@ -../../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/bfgs/water_64.xyz b/examples/temp/pes/geop/waterbox/bfgs/water_64.xyz new file mode 120000 index 000000000..683f52959 --- /dev/null +++ b/examples/temp/pes/geop/waterbox/bfgs/water_64.xyz @@ -0,0 +1 @@ +../../../../init_files/water_64.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/cg/input.xml b/examples/temp/pes/geop/waterbox/cg/input.xml index 11feb189b..e6ea9da5f 100644 --- a/examples/temp/pes/geop/waterbox/cg/input.xml +++ b/examples/temp/pes/geop/waterbox/cg/input.xml @@ -12,7 +12,7 @@ - water_216.xyz + water_64.xyz diff --git a/examples/temp/pes/geop/waterbox/cg/water_216.xyz b/examples/temp/pes/geop/waterbox/cg/water_216.xyz deleted file mode 120000 index 95291df24..000000000 --- a/examples/temp/pes/geop/waterbox/cg/water_216.xyz +++ /dev/null @@ -1 +0,0 @@ -../../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/cg/water_64.xyz b/examples/temp/pes/geop/waterbox/cg/water_64.xyz new file mode 120000 index 000000000..683f52959 --- /dev/null +++ b/examples/temp/pes/geop/waterbox/cg/water_64.xyz @@ -0,0 +1 @@ +../../../../init_files/water_64.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/lbfgs/input.xml b/examples/temp/pes/geop/waterbox/lbfgs/input.xml index 9d1eb692e..279dc0257 100644 --- a/examples/temp/pes/geop/waterbox/lbfgs/input.xml +++ b/examples/temp/pes/geop/waterbox/lbfgs/input.xml @@ -12,7 +12,7 @@ - water_216.xyz + water_64.xyz diff --git a/examples/temp/pes/geop/waterbox/lbfgs/water_216.xyz b/examples/temp/pes/geop/waterbox/lbfgs/water_216.xyz deleted file mode 120000 index 95291df24..000000000 --- a/examples/temp/pes/geop/waterbox/lbfgs/water_216.xyz +++ /dev/null @@ -1 +0,0 @@ -../../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/lbfgs/water_64.xyz b/examples/temp/pes/geop/waterbox/lbfgs/water_64.xyz new file mode 120000 index 000000000..683f52959 --- /dev/null +++ b/examples/temp/pes/geop/waterbox/lbfgs/water_64.xyz @@ -0,0 +1 @@ +../../../../init_files/water_64.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/sd/input.xml b/examples/temp/pes/geop/waterbox/sd/input.xml index 85c49914f..ff259cf1f 100644 --- a/examples/temp/pes/geop/waterbox/sd/input.xml +++ b/examples/temp/pes/geop/waterbox/sd/input.xml @@ -12,7 +12,7 @@ - water_216.xyz + water_64.xyz diff --git a/examples/temp/pes/geop/waterbox/sd/water_216.xyz b/examples/temp/pes/geop/waterbox/sd/water_216.xyz deleted file mode 120000 index 95291df24..000000000 --- a/examples/temp/pes/geop/waterbox/sd/water_216.xyz +++ /dev/null @@ -1 +0,0 @@ -../../../../init_files/water_216.xyz \ No newline at end of file diff --git a/examples/temp/pes/geop/waterbox/sd/water_64.xyz b/examples/temp/pes/geop/waterbox/sd/water_64.xyz new file mode 120000 index 000000000..683f52959 --- /dev/null +++ b/examples/temp/pes/geop/waterbox/sd/water_64.xyz @@ -0,0 +1 @@ +../../../../init_files/water_64.xyz \ No newline at end of file diff --git a/ipi/engine/atoms.py b/ipi/engine/atoms.py index 73b4c1ce4..8638a82c2 100644 --- a/ipi/engine/atoms.py +++ b/ipi/engine/atoms.py @@ -148,7 +148,7 @@ def __init__(self, natoms, _prebind=None): name="kstress", func=self.get_kstress, dependencies=[self._p, self._m] ) - def copy(self): + def clone(self): """Creates a new Atoms object. Returns: @@ -156,9 +156,9 @@ def copy(self): """ newat = Atoms(self.natoms) - newat.q[:] = self.q - newat.p[:] = self.p - newat.m[:] = self.m + newat.q[:] = dstrip(self.q) + newat.p[:] = dstrip(self.p) + newat.m[:] = dstrip(self.m) newat.names[:] = self.names return newat diff --git a/ipi/engine/cell.py b/ipi/engine/cell.py index ca69a20d2..e60eda0bb 100644 --- a/ipi/engine/cell.py +++ b/ipi/engine/cell.py @@ -14,7 +14,7 @@ from ipi.utils.mathtools import * -__all__ = ["Cell"] +__all__ = ["Cell", "GenericCell"] class Cell: @@ -95,15 +95,14 @@ def array_pbc(self, pos): system box. """ - s = dstrip(pos).copy() - s.shape = (len(pos) // 3, 3) + s = dstrip(pos).copy().reshape(-1, 3) s = np.dot(dstrip(self.ih), s.T) s = s - np.round(s) s = np.dot(dstrip(self.h), s).T - pos[:] = s.reshape((len(s) * 3)) + pos[:] = s.reshape(-1) def minimum_distance(self, atom1, atom2): """Takes two atoms and tries to find the smallest vector between two @@ -130,3 +129,49 @@ def minimum_distance(self, atom1, atom2): dproperties(Cell, ["h", "ih", "V"]) + + +class GenericCell(Cell): + """A cell class that does not assume upper-triangular values. + + Depend objects: + h: An array giving the lattice vector matrix. + ih: An array giving the inverse of the lattice vector matrix. + V: The volume of the cell. + """ + + def __init__(self, h=None): + """Initialises base cell class. + + Args: + h: Optional array giving the initial lattice vector matrix. The + reference cell matrix is set equal to this. + """ + + if h is None: + h = np.zeros((3, 3), float) + + self._h = depend_array(name="h", value=h) + self._ih = depend_array( + name="ih", + value=np.zeros((3, 3), float), + func=self.get_ih, + dependencies=[self._h], + ) + self._V = depend_value(name="V", func=self.get_volume, dependencies=[self._h]) + + def clone(self): + return GenericCell(dstrip(self.h).copy()) + + def get_ih(self): + """Inverts the lattice vector matrix.""" + + return np.linalg.inv(self.h) + + def get_volume(self): + """Calculates the volume of the system box.""" + + return np.linalg.det(self.h) + + +dproperties(GenericCell, ["h", "ih", "V"]) diff --git a/ipi/engine/forcefields.py b/ipi/engine/forcefields.py index d764902b6..bb3c026cd 100644 --- a/ipi/engine/forcefields.py +++ b/ipi/engine/forcefields.py @@ -13,9 +13,12 @@ import threading import json import sys +from contextlib import nullcontext import numpy as np +from ipi.engine.cell import GenericCell +from ipi.utils.prng import Random from ipi.utils.softexit import softexit from ipi.utils.messages import info, verbosity, warning from ipi.interfaces.sockets import InterfaceSocket @@ -24,6 +27,11 @@ from ipi.utils.units import unit_to_internal from ipi.utils.distance import vector_separation from ipi.pes import __drivers__ +from ipi.utils.mathtools import ( + get_rotation_quadrature_legendre, + get_rotation_quadrature_lebedev, + random_rotation, +) plumed = None @@ -67,7 +75,7 @@ class ForceField: def __init__( self, - latency=1.0, + latency=1e-4, offset=0.0, name="", pars=None, @@ -223,22 +231,26 @@ def _poll_loop(self): finished. """ - info(" @ForceField: Starting the polling thread main loop.", verbosity.low) + info( + f" @ForceField ({self.name}): Starting the polling thread main loop.", + verbosity.low, + ) while self._doloop[0]: time.sleep(self.latency) if len(self.requests) > 0: self.poll() - def release(self, request): + def release(self, request, lock=True): """Removes a request from the evaluation queue. Args: request: The id of the job to release. + lock: whether we should apply a threadlock here """ """Frees up a request.""" - with self._threadlock: + with self._threadlock if lock else nullcontext(): if request in self.requests: try: self.requests.remove(request) @@ -419,6 +431,8 @@ def __init__( super().__init__(latency, offset, name, pars, dopbc, active, threaded) + if pars is None: + pars = {} # defaults no pars self.pes = pes try: self.driver = __drivers__[self.pes](verbose=verbosity.high, **pars) @@ -1159,7 +1173,7 @@ def __init__( pars=pars, dopbc=dopbc, active=active, - threaded=True, + threaded=True, # hardcoded, otherwise won't work! ) if len(fflist) == 0: raise ValueError( @@ -1442,6 +1456,209 @@ def poll(self): r["status"] = "Done" +class FFRotations(ForceField): + """Forcefield to manipulate models that are not exactly rotationally equivariant. + Can be used to evaluate a different random rotation at each evaluation, or to average + over a regular grid of Euler angles""" + + def __init__( + self, + latency=1.0, + offset=0.0, + name="", + pars=None, + dopbc=True, + active=np.array([-1]), + threaded=True, + prng=None, + ffsocket=None, + ffdirect=None, + random=False, + inversion=False, + grid_order=1, + grid_mode="lebedev", + ): + super( + FFRotations, self + ).__init__( # force threaded execution to handle sub-ffield + latency, offset, name, pars, dopbc, active, threaded=True + ) + + if prng is None: + warning("No PRNG provided, will initialize one", verbosity.low) + self.prng = Random() + else: + self.prng = prng + self.ffsocket = ffsocket + self.ffdirect = ffdirect + if ffsocket is None or self.ffsocket.name == "__DUMMY__": + if ffdirect is None or self.ffdirect.name == "__DUMMY__": + raise ValueError( + "Must specify a non-default value for either `ffsocket` or `ffdirect` into `ffrotations`" + ) + else: + self.ff = self.ffdirect + elif ffdirect is None or self.ffdirect.name == "__DUMMY__": + self.ff = self.ffsocket + else: + raise ValueError( + "Cannot specify both `ffsocket` and `ffdirect` into `ffrotations`" + ) + + self.random = random + self.inversion = inversion + self.grid_order = grid_order + self.grid_mode = grid_mode + + if self.grid_mode == "lebedev": + self._rotations = get_rotation_quadrature_lebedev(self.grid_order) + elif self.grid_mode == "legendre": + self._rotations = get_rotation_quadrature_legendre(self.grid_order) + else: + raise ValueError(f"Invalid quadrature {self.grid_mode}") + info( + f""" +# Generating {self.grid_mode} rotation quadrature of order {self.grid_order}. +# Grid contains {len(self._rotations)} proper rotations. +{("# Inversion is also active, doubling the number of evaluations." if self.inversion else "")}""", + verbosity.low, + ) + + def bind(self, output_maker=None): + super().bind(output_maker) + self.ff.bind(output_maker) + + def start(self): + super().start() + self.ff.start() + + def queue(self, atoms, cell, reqid=-1): + + # launches requests for all of the rotations FF objects + ffh = [] # this is the list of "inner" FF requests + rots = [] # this is a list of tuples of (rotation matrix, weight) + if self.random: + R_random = random_rotation(self.prng, improper=True) + else: + R_random = np.eye(3) + + for R, w, _ in self._rotations: + R = R @ R_random + + rot_atoms = atoms.clone() + rot_cell = GenericCell( + R @ dstrip(cell.h).copy() + ) # NB we need generic cell orientation + rot_atoms.q[:] = (dstrip(rot_atoms.q).reshape(-1, 3) @ R.T).flatten() + + rots.append((R, w)) + ffh.append(self.ff.queue(rot_atoms, rot_cell, reqid)) + + if self.inversion: + # also add a "flipped rotation" to the evaluation list + R = R * -1 + + rot_cell = GenericCell(R @ dstrip(cell.h).copy()) + rot_atoms = atoms.clone() + rot_atoms.q[:] = (dstrip(rot_atoms.q).reshape(-1, 3) @ R.T).flatten() + + rots.append((R, w)) + ffh.append(self.ff.queue(rot_atoms, rot_cell, reqid)) + + # creates the request with the help of the base class, + # making sure it already contains a handle to the list of FF + # requests + req = ForceField.queue( + self, atoms, cell, reqid, template=dict(ff_handles=ffh, rots=rots) + ) + req["status"] = "Running" + req["t_dispatched"] = time.time() + return req + + def check_finish(self, r): + """Checks if all sub-requests associated with a given + request are finished""" + + for ff_r in r["ff_handles"]: + if ff_r["status"] != "Done": + return False + return True + + def gather(self, r): + """Collects results from all sub-requests, and assemble the committee of models.""" + + r["result"] = [ + 0.0, + np.zeros(len(r["pos"]), float), + np.zeros((3, 3), float), + "", + ] + + # list of pointers to the forcefield requests. shallow copy so we can remove stuff + rot_handles = r["ff_handles"].copy() + rots = r["rots"].copy() + + # Gathers the forcefield energetics and extras + pots = [] + frcs = [] + virs = [] + xtrs = [] + quad_w = [] + for ff_r, (R, w) in zip(rot_handles, rots): + pots.append(ff_r["result"][0]) + # must rotate forces and virial back into the original reference frame + frcs.append((ff_r["result"][1].reshape(-1, 3) @ R).flatten()) + virs.append((R.T @ ff_r["result"][2] @ R)) + xtrs.append(ff_r["result"][3]) + quad_w.append(w) + + quad_w = np.array(quad_w) + pots = np.array(pots) + frcs = np.array(frcs).reshape(len(frcs), -1) + virs = np.array(virs).reshape(-1, 3, 3) + + # Computes the mean energetics (using the quadrature weights) + mean_pot = np.sum(pots * quad_w, axis=0) / quad_w.sum() + mean_frc = np.sum(frcs * quad_w[:, np.newaxis], axis=0) / quad_w.sum() + mean_vir = ( + np.sum(virs * quad_w[:, np.newaxis, np.newaxis], axis=0) / quad_w.sum() + ) + + # Sets the output of the committee model. + r["result"][0] = mean_pot + r["result"][1] = mean_frc + r["result"][2] = mean_vir + r["result"][3] = { + "o3grid_pots": pots + } # this is the list of potentials on a grid, for monitoring + + # "dissolve" the extras dictionaries into a list + if isinstance(xtrs[0], dict): + for k in xtrs[0].keys(): + r["result"][3][k] = [] + for x in xtrs: + r["result"][3][k].append(x[k]) + else: + r["result"][3]["raw"] = [] + for x in xtrs: + r["result"][3]["raw"].append(x) + + for ff_r in r["ff_handles"]: + self.ff.release(ff_r) + + def poll(self): + """Polls the forcefield object to check if it has finished.""" + + with self._threadlock: + for r in self.requests: + if "ff_handles" in r and r["status"] != "Done" and self.check_finish(r): + r["t_finished"] = time.time() + self.gather(r) + r["result"][0] -= self.offset + r["status"] = "Done" + self.release(r, lock=False) + + class PhotonDriver: """ Photon driver for a single cavity mode diff --git a/ipi/inputs/forcefields.py b/ipi/inputs/forcefields.py index 5174a8c23..245908219 100644 --- a/ipi/inputs/forcefields.py +++ b/ipi/inputs/forcefields.py @@ -20,6 +20,7 @@ FFCommittee, FFdmd, FFCavPhSocket, + FFRotations, ) from ipi.interfaces.sockets import InterfaceSocket from ipi.pes import __drivers__ @@ -27,6 +28,8 @@ from ipi.inputs.initializer import * from ipi.utils.inputvalue import * from ipi.utils.messages import verbosity, warning +from ipi.utils.prng import Random +from ipi.inputs.prng import InputRandom __all__ = [ "InputFFSocket", @@ -39,6 +42,7 @@ "InputFFCommittee", "InputFFdmd", "InputFFCavPhSocket", + "InputFFRotations", ] @@ -263,7 +267,7 @@ def store(self, ff): ff: A ForceField object with a FFSocket forcemodel object. """ - if (not type(ff) is FFSocket) and (not type(ff) is FFCavPhSocket): + if type(ff) not in [FFSocket, FFCavPhSocket]: raise TypeError( "The type " + type(ff).__name__ + " is not a valid socket forcefield" ) @@ -471,7 +475,7 @@ class InputFFDebye(InputForceField): "default": input_default(factory=np.zeros, args=(0,)), "help": "Specifies the Hessian of the harmonic potential. " "Default units are atomic. Units can be specified only by xml attribute. " - r"Implemented options are: 'atomic_unit', 'ev/ang\^2'", + r"Implemented options are: 'atomic_unit', 'ev/ang^2'", "dimension": "hessian", }, ), @@ -921,6 +925,127 @@ def fetch(self): ) +class InputFFRotations(InputForceField): + default_help = """ +Wraps around another forcefield to evaluate it over one or more +rotated copies of the physical system. This is useful when +interacting with models that are not exactly invariant/covariant +with respect to rigid rotations. +Besides the parameters defining how averaging is to be performed +(using an integration grid, and/or randomizing the orientation at +each step) the should contain either a +or a block that computes the "base" model. Note that +this forcefield should be given a separate `name`, but that you +cannot access this "inner" forcefield from other parts of the +input file. +""" + default_label = "FFROTATIONS" + + fields = copy(InputForceField.fields) + + fields["random"] = ( + InputValue, + { + "dtype": bool, + "default": False, + "help": """Applies a random rotation at each evaluation. """, + }, + ) + + fields["grid_order"] = ( + InputValue, + { + "dtype": int, + "default": 1, + "help": """Sums over a grid of rotations of the given order. + Note that the number of rotations increases rapidly with the order: + e.g. for a legendre grid + '1' leads to a single rotation, '2' to 18, '3' to 75 rotations, while + for a lebedev grid '3' contains 18 rotations, '5' 70 rotations and so on. + """, + }, + ) + + fields["grid_mode"] = ( + InputValue, + { + "dtype": str, + "options": ["legendre", "lebedev"], + "default": "legendre", + "help": """Defines the type of integration grid. + Lebedev grids are usually more efficient in integrating. + """, + }, + ) + + fields["inversion"] = ( + InputValue, + { + "dtype": bool, + "default": False, + "help": """Always applies the improper version of each rotation in the + grid (or the randomly-sampled rotation). Doubles the evaluations and makes + the model exactly equivariant to inversion.""", + }, + ) + + fields["prng"] = ( + InputRandom, + { + "help": InputRandom.default_help, + "default": input_default(factory=Random), + }, + ) + + fields["ffsocket"] = ( + InputFFSocket, + { + "help": InputFFSocket.default_help, + "default": input_default(factory=FFSocket, kwargs={"name": "__DUMMY__"}), + }, + ) + + fields["ffdirect"] = ( + InputFFDirect, + { + "help": InputFFDirect.default_help, + "default": input_default(factory=FFDirect, kwargs={"name": "__DUMMY__"}), + }, + ) + + def store(self, ff): + """Store the base and rotation quadrature parameters""" + + super(InputFFRotations, self).store(ff) + self.inversion.store(ff.inversion) + self.grid_order.store(ff.grid_order) + self.grid_mode.store(ff.grid_mode) + self.random.store(ff.random) + self.prng.store(ff.prng) + self.ffsocket.store(ff.ffsocket) + self.ffdirect.store(ff.ffdirect) + + def fetch(self): + """Fetches all of the FF objects""" + + return FFRotations( + pars=self.parameters.fetch(), + name=self.name.fetch(), + latency=self.latency.fetch(), + offset=self.offset.fetch(), + dopbc=self.pbc.fetch(), + active=self.activelist.fetch(), + threaded=self.threaded.fetch(), + prng=self.prng.fetch(), + ffsocket=self.ffsocket.fetch(), + ffdirect=self.ffdirect.fetch(), + grid_order=self.grid_order.fetch(), + grid_mode=self.grid_mode.fetch(), + random=self.random.fetch(), + inversion=self.inversion.fetch(), + ) + + class InputFFCavPhSocket(InputFFSocket): default_help = """A cavity molecular dynamics driver for vibraitonal strong coupling. In the current implementation, only a single cavity mode polarized along the x and y directions is coupled to the molecules. diff --git a/ipi/inputs/simulation.py b/ipi/inputs/simulation.py index 796316d59..3c1c4054f 100644 --- a/ipi/inputs/simulation.py +++ b/ipi/inputs/simulation.py @@ -187,6 +187,10 @@ class InputSimulation(Input): iforcefields.InputFFCommittee, {"help": iforcefields.InputFFCommittee.default_help}, ), + "ffrotations": ( + iforcefields.InputFFRotations, + {"help": iforcefields.InputFFRotations.default_help}, + ), "ffsgdml": ( iforcefields.InputFFsGDML, {"help": iforcefields.InputFFsGDML.default_help}, @@ -246,7 +250,7 @@ def store(self, simul): _obj, ) in enumerate(_fflist + simul.syslist): if self.extra[_ii] == 0: - if isinstance(_obj, eforcefields.FFSocket): + if type(_obj) is eforcefields.FFSocket: _iobj = iforcefields.InputFFSocket() _iobj.store(_obj) self.extra[_ii] = ("ffsocket", _iobj) @@ -282,6 +286,10 @@ def store(self, simul): _iobj = iforcefields.InputFFCommittee() _iobj.store(_obj) self.extra[_ii] = ("ffcommittee", _iobj) + elif isinstance(_obj, eforcefields.FFRotations): + _iobj = iforcefields.InputFFRotations() + _iobj.store(_obj) + self.extra[_ii] = ("ffrotations", _iobj) elif isinstance(_obj, eforcefields.FFCavPhSocket): _iobj = iforcefields.InputFFCavPhSocket() _iobj.store(_obj) @@ -340,6 +348,7 @@ def fetch(self): "ffsgdml", "ffyaff", "ffcommittee", + "ffrotations", "ffcavphsocket", ]: new_ff = v.fetch() diff --git a/ipi/utils/__init__.py b/ipi/utils/__init__.py index 24b383cd3..808781bfd 100644 --- a/ipi/utils/__init__.py +++ b/ipi/utils/__init__.py @@ -16,6 +16,7 @@ "messages", "softexit", "setup", + "tools", "io", "constrtools", ] diff --git a/ipi/utils/lebedev_grids.npy b/ipi/utils/lebedev_grids.npy new file mode 100644 index 000000000..e8c36c3be Binary files /dev/null and b/ipi/utils/lebedev_grids.npy differ diff --git a/ipi/utils/mathtools.py b/ipi/utils/mathtools.py index 000043233..e6a6341b5 100644 --- a/ipi/utils/mathtools.py +++ b/ipi/utils/mathtools.py @@ -6,6 +6,7 @@ import math +from importlib import resources import numpy as np @@ -25,6 +26,8 @@ "logsumlog", "sinch", "gaussian_inv", + "get_rotation_quadrature_legendre", + "get_rotation_quadrature_lebedev", ] @@ -521,3 +524,127 @@ def f(w): )[0] return 2 * integral / np.pi + + +def quat2rot(q): + """ + Convert a normalized quaternion into a 3x3 rotation matrix. + + Args: + q : list or array-like + A normalized quaternion [q1, q2, q3, q4] where q1 is the scalar part. + + Returns: + numpy.ndarray + The corresponding 3x3 rotation matrix. + """ + + q1, q2, q3, q4 = q + + # Compute the rotation matrix elements + R = np.array( + [ + [1 - 2 * (q3**2 + q4**2), 2 * (q2 * q3 - q1 * q4), 2 * (q2 * q4 + q1 * q3)], + [2 * (q2 * q3 + q1 * q4), 1 - 2 * (q2**2 + q4**2), 2 * (q3 * q4 - q1 * q2)], + [2 * (q2 * q4 - q1 * q3), 2 * (q3 * q4 + q1 * q2), 1 - 2 * (q2**2 + q3**2)], + ] + ) + + return R + + +def euler_zxz_to_matrix(theta, v, w): + """ + Generate a rotation matrix from Euler angles [theta, v, w] in the ZXZ convention. + + Args: + theta, v, w : float + Euler angles in radians. + + Returns: + numpy.ndarray + The corresponding 3x3 rotation matrix. + """ + + # Rotation matrix around the Z-axis by theta + R_z_theta = np.array( + [ + [np.cos(theta), -np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 0], + [0, 0, 1], + ] + ) + + # Rotation matrix around the X-axis by v + R_x_v = np.array([[1, 0, 0], [0, np.cos(v), -np.sin(v)], [0, np.sin(v), np.cos(v)]]) + + # Rotation matrix around the Z-axis by w + R_z_w = np.array([[np.cos(w), -np.sin(w), 0], [np.sin(w), np.cos(w), 0], [0, 0, 1]]) + + # Total rotation matrix is the product of these matrices: R_z(w) * R_x(v) * R_z(theta) + rotation_matrix = R_z_w @ R_x_v @ R_z_theta + + return rotation_matrix + + +def roots_legendre(L): + """Replicates scipy.special.roots_legendre using only numpy""" + + legendre_poly = np.polynomial.legendre.Legendre.basis(L) + roots = np.polynomial.legendre.legroots(legendre_poly.coef) + legendre_poly_deriv = legendre_poly.deriv() + + # Calculate weights using the formula + weights = 2 / ((1 - roots**2) * (legendre_poly_deriv(roots) ** 2)) + + return roots, weights + + +def get_rotation_quadrature_legendre(L): + if L == 1: + # returns the identity (for some reason this algo below generates a different rotation) + return [(np.eye(3), 2.0, [0, 0, 0])] + quads = [] + for theta_index in range(0, 2 * L - 1): + for w_index in range(0, 2 * L - 1): + theta = 2 * np.pi * theta_index / (2 * L - 1) + w = 2 * np.pi * w_index / (2 * L - 1) + roots_legendre_now, weights_now = roots_legendre(L) + all_v = np.arccos(roots_legendre_now) + for v, weight in zip(all_v, weights_now): + angles = [theta, v, w] + rotation_matrix = euler_zxz_to_matrix(*angles) + quads.append((rotation_matrix, weight, angles)) + + return quads + + +def get_rotation_quadrature_lebedev(L): + with resources.path("ipi.utils", "lebedev_grids.npy") as file_path: + lebedev = np.load(file_path, allow_pickle=True).item() + if not L in lebedev: + raise ValueError(f"There is no Lebedev grid of order {L} available") + leb_quad = lebedev[L] + quads = [] + for theta_index in range(0, L): + theta = 2 * np.pi * theta_index / L + for w, v, weight in leb_quad: + angles = [w, v, theta] + rotation_matrix = euler_zxz_to_matrix(*angles) + quads.append((rotation_matrix, weight, angles)) + return quads + + +def random_rotation(prng, improper=True): + """Generates a (uniform) random rotation matrix""" + + quaternion = prng.gvec(shape=(4,)) + quaternion /= np.sqrt((quaternion**2).sum()) + + rot = quat2rot(quaternion) + + # randomly generate an improper rotation + if improper and prng.u < 0.5: + rot *= -1.0 + + return rot diff --git a/ipi/utils/setup.py b/ipi/utils/setup.py index d8fd38ab2..264ec43ce 100644 --- a/ipi/utils/setup.py +++ b/ipi/utils/setup.py @@ -64,7 +64,13 @@ def install_driver(force_install=False): try: subprocess.run(["git", "init"], cwd=temp_dir) subprocess.run( - ["git", "remote", "add", "origin", "https://github.com/i-pi/i-pi.git"], + [ + "git", + "remote", + "add", + "origin", + "https://github.com/lab-cosmo/i-pi.git", + ], cwd=temp_dir, ) subprocess.run( @@ -74,7 +80,9 @@ def install_driver(force_install=False): os.path.join(temp_dir, ".git", "info", "sparse-checkout"), "w" ) as f: f.write("drivers/f90\n") - subprocess.run(["git", "pull", "--depth=1", "origin", "main"], cwd=temp_dir) + subprocess.run( + ["git", "pull", "--depth=1", "origin", "feat/equivator"], cwd=temp_dir + ) except: warning( "Failed to fetch the drivers folder from i-PI github repository", diff --git a/ipi_tests/examples/exampletools.py b/ipi_tests/examples/exampletools.py index a78ac3c5d..f29ef10ae 100644 --- a/ipi_tests/examples/exampletools.py +++ b/ipi_tests/examples/exampletools.py @@ -55,5 +55,6 @@ def create_client_list(self, driver_info, nid, test_settings): test_settings, ) return clients - except: - raise RuntimeError("Couldn't modify the xml file") + except Exception as ex: + print("Couldn't modify the xml file") + raise ex diff --git a/ipi_tests/test_tools.py b/ipi_tests/test_tools.py index b4f4dff35..f43b25ecf 100644 --- a/ipi_tests/test_tools.py +++ b/ipi_tests/test_tools.py @@ -9,6 +9,7 @@ clean_all = False debug = False +TIMEOUT = 20 fortran_driver_models = [ "dummy", @@ -31,6 +32,7 @@ "doublewell_1D", "doublewell", "gas", + "noo3-h2o", ] # We should do this automatically but for now we do it explicitly here @@ -185,10 +187,11 @@ def modify_xml_4_dummy_test( root = tree.getroot() clients = list() + ff_roots = [root] if len(root.findall("ffcommittee")) > 0: - ff_roots = root.findall("ffcommittee") - else: - ff_roots = [root] + ff_roots += root.findall("ffcommittee") + if len(root.findall("ffrotations")) > 0: + ff_roots += root.findall("ffrotations") ff_sockets = [] for ff_root in ff_roots: @@ -322,7 +325,7 @@ def run(self, cwd, nid): print("List all files /tmp/ipi_*") for filename in glob.glob("/tmp/ipi_*"): print(filename) - ipi_error = ipi.communicate(timeout=120)[1].decode("ascii") + ipi_error = ipi.communicate(timeout=TIMEOUT)[1].decode("ascii") print(ipi_error) return "Could not find the i-PI UNIX socket" @@ -375,13 +378,13 @@ def run(self, cwd, nid): drivers.append(driver) # check i-pi errors - ipi_out, ipi_error = ipi.communicate(timeout=60) + ipi_out, ipi_error = ipi.communicate(timeout=TIMEOUT) assert ipi.returncode == 0, "i-PI error occurred: {}".format(ipi_error) # check driver errors for driver in drivers: # if i-PI has ended, we can wait for the driver to quit - driver_out, driver_err = driver.communicate(timeout=60) + driver_out, driver_err = driver.communicate(timeout=TIMEOUT) assert ( driver.returncode == 0 ), "Driver error occurred: {}\n Driver Output: {}".format( diff --git a/setup.cfg b/setup.cfg index eda5e7f97..bf28c22a5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,10 +21,10 @@ classifiers = Operating System :: POSIX :: Linux License :: OSI Approved :: GNU General Public License v3 (GPLv3) Programming Language :: Python :: 3 - Programming Language :: Python :: 3.6 Programming Language :: Python :: 3.7 Programming Language :: Python :: 3.8 Programming Language :: Python :: 3.9 + Programming Language :: Python :: 3.10 Topic :: Scientific/Engineering :: Chemistry Topic :: Scientific/Engineering :: Physics diff --git a/setup.py b/setup.py index e865e9c79..a4d3c095f 100644 --- a/setup.py +++ b/setup.py @@ -7,5 +7,8 @@ setup( packages=[*find_packages(exclude=["ipi_tests*"])], + package_data={ + "ipi.utils": ["lebedev_grids.npy"], # Include the npy file here + }, scripts=scripts, )