Skip to content

Commit

Permalink
Merge branch 'master' into driver-auto-pes
Browse files Browse the repository at this point in the history
  • Loading branch information
mahrossi authored May 8, 2024
2 parents fc85e2c + 13943fa commit 1431339
Show file tree
Hide file tree
Showing 267 changed files with 16,944 additions and 14,107 deletions.
7 changes: 5 additions & 2 deletions docs/src/features.rst
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,12 @@ Fermionic statistics can be obtained by a reweigthing procedure by post
processing the simulation (see ref. 2 below).

| **Main contributors:** Yotam Feldman, Barak Hirshberg
| **Implementation:**
| **Theory and Implementation:**
| Y.M.Y. Feldman and B. Hirshberg *“Quadratic scaling bosonic path
integral molecular dynamics simulations”*, in preparation.
integral molecular dynamics simulations”*, J. Chem. Phys. 159, 154107 (2023) DOI:
`10.1063/5.0173749 <https://doi.org/10.1063/5.0173749>`__
— BibTeX:
`fetch <https://pubs.aip.org/Citation/Download?resourceId=2917443&resourceType=3&citationFormat=2>`__
| **Theory:**
| B. Hirshberg, V. Rizzi and M. Parrinello, *“Path integral molecular
dynamics for bosons”*, Proc. Natl. Acad. Sci. U.S.A. 116 (43)
Expand Down
82 changes: 82 additions & 0 deletions drivers/f90/LJ.f90
Original file line number Diff line number Diff line change
Expand Up @@ -209,4 +209,86 @@ SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list,

END SUBROUTINE

SUBROUTINE LJMix_getall(n_type2, rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
! Calculates the LJ potential energy and virial and the forces
! acting on all the atoms, using 0.4 sigma and 0.4 epsilon for the first n_type2 atoms
!
! Args:
! n_type2: First atom are considered as type 2
! rc: The cut-off radius.
! sigma: The LJ distance parameter. NOTE: x0 of the spring is set to be equal to sigma
! eps: The LJ energy parameter. NOTE: stiffness of the spring is set to be equal: 36*2^(2/3) eps
! natoms: The number of atoms in the system.
! atoms: A vector holding all the atom positions.
! cell_h: The simulation box cell vector matrix.
! cell_ih: The inverse of the simulation box cell vector matrix.
! index_list: A array giving the last index of n_list that
! gives the neighbours of a given atom.
! n_list: An array giving the indices of the atoms that neighbour
! the atom determined by index_list.
! pot: The total potential energy of the system.
! forces: An array giving the forces acting on all the atoms.
! virial: The virial tensor, not divided by the volume.

DOUBLE PRECISION, INTENT(IN) :: rc
DOUBLE PRECISION, INTENT(IN) :: sigma
DOUBLE PRECISION, INTENT(IN) :: eps
INTEGER, INTENT(IN) :: natoms
DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(IN) :: atoms
DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h
DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih
INTEGER, DIMENSION(natoms), INTENT(IN) :: index_list
INTEGER, DIMENSION(natoms*(natoms-1)/2), INTENT(IN) :: n_list
DOUBLE PRECISION, INTENT(OUT) :: pot
DOUBLE PRECISION, DIMENSION(natoms,3), INTENT(OUT) :: forces
DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: virial

INTEGER i, j, k, l, start, n_type2
DOUBLE PRECISION, DIMENSION(3) :: fij, rij
DOUBLE PRECISION sigmaij, epsij, r2, pot_ij, pot_lr, vir_lr, volume

start = 1
DO i = 1, natoms - 1
! Only loops over the neigbour list, not all the atoms.
DO j = start, index_list(i)
CALL vector_separation(cell_h, cell_ih, atoms(i,:), atoms(n_list(j),:), rij, r2)
IF (r2 < rc*rc) THEN ! Only calculates contributions between neighbouring particles.
sigmaij = 1 ! geometric average of interaction parameters
epsij = 1
IF (i .le. n_type2) THEN
sigmaij = 0.6
epsij = 0.6
ENDIF
IF (i .le. n_type2) THEN
sigmaij = sigmaij*0.6
epsij = epsij*0.6
ENDIF
sigmaij = sigma*DSQRT(sigmaij)
epsij = eps*DSQRT(epsij)

CALL LJ_fij(sigmaij, epsij, rij, sqrt(r2), pot_ij, fij)

forces(i,:) = forces(i,:) + fij
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.
virial(k,l) = virial(k,l) + fij(k)*rij(l)
ENDDO
ENDDO
ENDIF
ENDDO
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)
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
virial(k,k) = virial(k,k) + vir_lr
ENDDO

END SUBROUTINE
END MODULE
15 changes: 9 additions & 6 deletions drivers/f90/driver.f90
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,8 @@ PROGRAM DRIVER
vstyle = 29
ELSEIF (trim(cmdbuffer) == "water_dip_pol") THEN
vstyle = 31
ELSEIF (trim(cmdbuffer) == "ljmix") THEN
vstyle = 32
ELSEIF (trim(cmdbuffer) == "qtip4pf-c-1") THEN
vstyle = 60
ELSEIF (trim(cmdbuffer) == "qtip4pf-c-2") THEN
Expand All @@ -193,7 +195,7 @@ PROGRAM DRIVER
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|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] "
STOP "ENDED"
ENDIF
ELSEIF (ccmd == 4) THEN
Expand Down Expand Up @@ -294,7 +296,6 @@ PROGRAM DRIVER
STOP "ENDED"
ENDIF
isinit = .true.

ELSEIF (23 == vstyle) THEN !MB
IF (par_count == 0) THEN ! defaults values
vpars(1) = 0.004737803248674678
Expand Down Expand Up @@ -337,10 +338,10 @@ PROGRAM DRIVER
STOP "ENDED"
ENDIF
isinit = .true.
ELSEIF (22 == vstyle) THEN !ljpolymer
ELSEIF (22 == vstyle .or. 32 == vstyle) THEN !ljpolymer or ljmix
IF (4/= par_count) THEN
WRITE(*,*) "Error: parameters not initialized correctly."
WRITE(*,*) "For ljpolymer potential use n_monomer,sigma,epsilon,cutoff"
WRITE(*,*) "For ljpolymer and ljmix potential use n_monomer,sigma,epsilon,cutoff"
STOP "ENDED"
ELSE
n_monomer = nint(vpars(1))
Expand Down Expand Up @@ -901,6 +902,8 @@ PROGRAM DRIVER
CALL SG_getall(rc, nat, atoms, cell_h, cell_ih, index_list, n_list, pot, forces, virial)
ELSEIF (vstyle == 22) THEN ! ljpolymer potential.
CALL ljpolymer_getall(n_monomer,rc,sigma,eps,stiffness,nat,atoms,cell_h,cell_ih,index_list,n_list,pot,forces,virial)
ELSEIF (vstyle == 32) THEN ! lj mixture.
CALL ljmix_getall(n_monomer,rc,sigma,eps,nat,atoms,cell_h,cell_ih,index_list,n_list,pot,forces,virial)
ENDIF
IF (verbose > 0) WRITE(*,*) " Calculated energy is ", pot
ENDIF
Expand Down Expand Up @@ -1048,15 +1051,15 @@ PROGRAM DRIVER
CONTAINS
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|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]"
WRITE(*,*) " -o 'comma_separated_parameters' [-v] "
WRITE(*,*) ""
WRITE(*,*) " For LJ potential use -o sigma,epsilon,cutoff "
WRITE(*,*) " For SG potential use -o cutoff "
WRITE(*,*) " For 1D/3D harmonic oscillator use -o k "
WRITE(*,*) " For 1D morse oscillators use -o r0,D,a"
WRITE(*,*) " For qtip4pf-efield use -o Ex,Ey,Ez with Ei in V/nm"
WRITE(*,*) " For ljpolymer use -o n_monomer,sigma,epsilon,cutoff "
WRITE(*,*) " For ljpolymer or lkmix use -o n_monomer,sigma,epsilon,cutoff "
WRITE(*,*) " For gas, dummy, use the optional -o sleep_seconds to add a delay"
WRITE(*,*) " For the ideal, qtip4pf*, zundel, ch4hcbe, nasa, doublewell or doublewell_1D no options are needed! "
END SUBROUTINE helpmessage
Expand Down
10 changes: 8 additions & 2 deletions drivers/f90/sockets.c
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ Can be linked to a FORTRAN code that does not support sockets natively.
#include <sys/types.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <netinet/tcp.h>
#include <sys/un.h>
#include <netdb.h>

Expand All @@ -67,7 +68,7 @@ ignored here for C compatibility.
*/

{
int sockfd, ai_err;
int sockfd, ai_err, flag;

if (*inet>0)
{ // creates an internet socket
Expand All @@ -88,7 +89,12 @@ ignored here for C compatibility.
// creates socket
sockfd = socket(res->ai_family, res->ai_socktype, res->ai_protocol);
if (sockfd < 0) { perror("Error opening socket"); exit(-1); }


// set TCP_NODELAY=1 to disable Nagle's algorithm as it slows down the small transactions for i-PI
flag = 1;
int result = setsockopt(sockfd, IPPROTO_TCP, TCP_NODELAY, (char *)&flag, sizeof(int));
if (result < 0) { perror("Error setting TCP_NODELAY"); }

// makes connection
if (connect(sockfd, res->ai_addr, res->ai_addrlen) < 0)
{ perror("Error opening INET socket: wrong port or server unreachable"); exit(-1); }
Expand Down
2 changes: 2 additions & 0 deletions drivers/py/driver.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,8 @@ def run_driver(
sock.connect("/tmp/ipi_" + address)
else:
sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
# this reduces latency for the small messages passed by the i-PI protocol
sock.setsockopt(socket.IPPROTO_TCP, socket.TCP_NODELAY, 1)
sock.connect((address, port))

f_init = False
Expand Down
32 changes: 32 additions & 0 deletions drivers/py/pes/elphmod.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
"""Interface with [elphmod](https://github.com/janberges/elphmod) MD driver."""

import sys
from .dummy import Dummy_driver

__DRIVER_NAME__ = "elphmod"
__DRIVER_CLASS__ = "ModelIIIDriver"

ERROR_MSG = """
A pickled driver instance is required (see elphmod.md.Driver.save).
Example: python3 driver.py -u -m elphmod -o driver.pickle
"""


class ModelIIIDriver(Dummy_driver):
"""Wrapper around elphmod MD driver."""

def check_arguments(self):
"""Check arguments and load driver instance."""

import elphmod

if len(self.args) != 1:
sys.exit(ERROR_MSG)

self.driver = elphmod.md.Driver.load(self.args[0])

def __call__(self, cell, pos):
"""Calculate energy and forces for given structure."""

return self.driver(cell, pos)
10 changes: 10 additions & 0 deletions examples/clients/elphmod/.gitignore
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
*.Wmat
*.epmatwp
*.ifc
*.pdf
*.pickle
*.wigner
*.xyz
*_hr.dat
RESTART
!driver.xyz
38 changes: 38 additions & 0 deletions examples/clients/elphmod/Makefile
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
# Makefile for the elphmod example
#
# This file is part of i-PI.
# i-PI Copyright (C) 2014-2024 i-PI developers
# See the "licenses" directory for full license information.

.PHONY: all model driver run run2 plot clean
all: plot

IPI = i-pi
MODEL = model_hr.dat model.ifc model.epmatwp model.wigner
DRIVER = driver.pickle driver.xyz
RUN = run.out run.pos_0.xyz run.for_0.xyz
PLOT = plot.pdf

model $(MODEL): model.py
python3 $<

driver $(DRIVER): driver.py $(MODEL)
python3 $<

run $(RUN): input.xml $(DRIVER)
$(IPI) $< &
sleep 5
i-pi-driver-py -u -m elphmod -o $(word 2,$^)
wait

run2: input.xml run.py $(DRIVER)
$(IPI) $< &
sleep 5
python3 $(word 2,$^)
wait

plot $(PLOT): plot.py $(RUN)
python3 $<

clean:
rm -rf $(MODEL) $(DRIVER) $(RUN) $(PLOT) RESTART EXIT '#'*
51 changes: 51 additions & 0 deletions examples/clients/elphmod/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
elphmod - i-PI example
======================

This example connects i-PI with [elphmod](https://github.com/janberges/elphmod),
which provides energy surfaces for distorted systems based on *ab initio* data
for the symmetric one, using the approximation of a linearized electron-lattice
coupling. This realizes *model III* by Schobert *et al.*, [SciPost Phys. **16**,
046 (2024)](https://doi.org/10.21468/SciPostPhys.16.2.046).

The Python packages `elphmod` and `matplotlib` (for plotting only) are required.
Both can be installed with `pip` or from `conda-forge`.

This example performs the structural optimization of a charge-density wave in
monolayer tantalum disulfide. Running `make` performs the following steps:

1. Usually, the calculation is based on a Wannier Hamiltonian from Wannier90
(`model_hr.dat`), force constants from the PHonon code of Quantum ESPRESSO
(`model.ifc`), and the corresponding representation of the coupling from a
[patched version](https://github.com/janberges/elphmod/tree/master/patches)
of the EPW code (`model.epmatwp`, `model.wigner`). To avoid the *ab initio*
step, here we use a nearest-neighbor model instead, which can be created with
`make model` or `python3 model.py`.

2. Based on these data, we can build a supercell model and set up the driver
with `make driver` or `python3 driver.py`. This will create files with the
driver object (`driver.pickle`) and initial ionic positions (`driver.xyz`).
Note that the dummy `driver.xyz` required for testing will be overwritten.
Git will show the file as modified (or deleted after running `make clean`).

3. Now we are ready to run i-PI and the driver, which will communicate via a
Unix domain socket, and perform the optimization with `make run`. This will
start i-PI in the background, wait until it is ready, and start the driver
with `i-pi-driver-py -u -m elphmod -o driver.pickle`. i-PI will output the
potential energies (`run.out`), positions (`run.pos_0.xyz`), and forces
(`run.for_0.xyz`) for all steps.

4. Finally, we can view the trajectory with `make plot` or `python3 plot.py`.
This should open an interactive figure and finally create `plot.pdf`, where
arrows indicate atomic displacements. Since we have used a simplistic model,
the periodic lattice distortion differs from the expected 3-by-3 pattern.

Note that saving the driver to an intermediate file `driver.pickle` has two
advantages: First, when running it through `i-pi-driver-py`, only a single
parameter (the filename) has to be specified on the command line, rather than
all details of the flexible setup in `driver.py`. Second, setting up the driver
can be quite time consuming, so it is more efficient to do this only once when
launching it multiple times later.

It is alternatively possible to run the driver through the module `ipi._driver`,
which is however only available if i-PI has been installed using `setuptools`,
e.g., through `pip`. This can be tested with `make run2` or `python3 run.py`.
26 changes: 26 additions & 0 deletions examples/clients/elphmod/driver.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
#!/usr/bin/env python3

"""See [documentation](https://janberges.github.io/elphmod/modules/md.html)."""

import elphmod

el = elphmod.el.Model("model", rydberg=True)
ph = elphmod.ph.Model("model.ifc", divide_mass=False)
elph = elphmod.elph.Model("model.epmatwp", "model.wigner", el, ph, divide_mass=False)

driver = elphmod.md.Driver(
elph,
nk=(12, 12),
nq=(2, 2),
supercell=(9, 9),
kT=0.02,
f=elphmod.occupations.marzari_vanderbilt,
n=1.0,
)

driver.kT = 0.005
driver.f = elphmod.occupations.fermi_dirac
driver.random_displacements()

driver.save("driver.pickle")
driver.to_xyz("driver.xyz")
5 changes: 5 additions & 0 deletions examples/clients/elphmod/driver.xyz
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
3
# CELL{H}: 6.3116853 -3.15584265 0 0 5.466079812 0 0 0 28.34589207
Ta 0.000000000 3.644053207 0.000000000
S 3.155842650 1.822026604 2.834589207
S 3.155842650 1.822026604 -2.834589207
22 changes: 22 additions & 0 deletions examples/clients/elphmod/input.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
<simulation>
<output prefix='run'>
<properties stride='1' filename='out'>[step, potential]</properties>
<trajectory stride='1' filename='pos'>positions</trajectory>
<trajectory stride='1' filename='for'>forces</trajectory>
</output>
<ffsocket name='elphmod' mode='unix' pbc='false'>
<address>localhost</address>
</ffsocket>
<system>
<initialize nbeads='1'>
<file mode='xyz'>driver.xyz</file>
</initialize>
<forces>
<force forcefield='elphmod'></force>
</forces>
<motion mode='minimize'>
<optimizer mode='bfgs'></optimizer>
<fixcom>True</fixcom>
</motion>
</system>
</simulation>
Loading

0 comments on commit 1431339

Please sign in to comment.