Merge branch 'master' into bug_fix_instanton

.gitignore

@@ -3,6 +3,7 @@
 bin/i-pi-driver
 doc/manual.pdf
 doc/manual.xml
+doc/latex/input_ref_idx.tex
 *.pyc
 *.pyo
 *.x

docs/src/features.rst

@@ -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)

drivers/f90/LJ.f90

@@ -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

drivers/f90/driver.f90

@@ -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
@@ -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
@@ -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
@@ -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))
@@ -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
@@ -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

drivers/f90/sockets.c

@@ -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>
 
@@ -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
@@ -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); }

drivers/py/driver.py

@@ -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