Merge branch 'master' into feat/socket_prefix

.github/workflows/lint.yml

@@ -18,5 +18,6 @@ jobs:
         uses: rbialon/flake8-annotations@v1
       - name: flake8 for syntax
         run: |
-           python -m pip install flake8==6.1.0  
-           flake8 --select=F --ignore= --ignore=F403,F405 --per-file-ignores=**__init__.py:F401 --statistics
+           python -m pip install flake8
+           flake8 --version
+           flake8 --select=F --ignore=F403,F405 --per-file-ignores=**__init__.py:F401 --statistics

.gitignore

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

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

@@ -177,6 +177,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
@@ -195,7 +197,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
@@ -298,7 +300,6 @@ PROGRAM DRIVER
             STOP "ENDED"
          ENDIF
          isinit = .true.
-
       ELSEIF (23 == vstyle) THEN !MB
          IF (par_count == 0) THEN ! defaults values
             vpars(1) = 0.004737803248674678
@@ -341,10 +342,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))
@@ -905,6 +906,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
@@ -1052,15 +1055,16 @@ 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' [-S sockets_prefix] [-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/py/pes/doubledoublewell.py

@@ -14,7 +14,6 @@
 __DRIVER_NAME__ = "double_double_well"
 __DRIVER_CLASS__ = "DDW_with_explicit_bath_driver"
 
-# np.set_printoptions(precision=14, suppress=True,threshold='nan',linewidth=1000)
 
 invcm2au = units.unit_to_internal("frequency", "inversecm", 1.0)
 A2au = units.unit_to_internal("length", "angstrom", 1.0)
@@ -54,16 +53,16 @@ def __init__(self, args=None, verbose=None):
         self.error_msg = """\nDW+explicit_bath driver expects 11 arguments.\n
         Example: python driver.py -m DoubleWell_with_explicit_bath -o wb1 (cm^-1) V1 (cm^-1) wb2 (cm^-1) V2 (cm^-1) coupling(au) mass delta(\AA) eta0  eps1 eps2  deltaQ omega_c(cm^-1)     \n
         python driver.py -m DoubleWell -o 500,2085,500,2085,0.1,1837,0.00,1,0,0,1,500\n"""
-        self.args = args
-        self.check_arguments()
+        super(DDW_with_explicit_bath_driver, self).__init__(
+            args, error_msg=self.error_msg
+        )
         self.init = False
 
     def check_arguments(self):
         """Function that checks the arguments required to run the driver"""
 
         try:
-            arglist = self.args.split(",")
-            param = list(map(float, arglist))
+            param = list(map(float, self.args))
             assert len(param) == 12
             wb1 = param[0] * invcm2au
             v1 = param[1] * invcm2au
@@ -109,8 +108,6 @@ def __call__(self, cell, pos):
         y = pos.reshape(-1, 1)[self.ndof // 2 + 1 :]
         assert self.ndof == q1.size + q2.size + x.size + y.size
 
-        # fq1 = np.zeros(q1.shape)
-        # fq2 = np.zeros(q2.shape)
         fx = np.zeros(x.shape)
         fy = np.zeros(x.shape)
 

drivers/py/pes/doublewell.py

@@ -15,8 +15,6 @@
 __DRIVER_CLASS__ = "DoubleWell_driver"
 
 
-# np.set_printoptions(precision=14, suppress=True,threshold='nan',linewidth=1000)
-
 invcm2au = units.unit_to_internal("frequency", "inversecm", 1.0)
 A2au = units.unit_to_internal("length", "angstrom", 1.0)
 
@@ -33,7 +31,7 @@ class DoubleWell_driver(Dummy_driver):
     def __init__(self, args=None, verbose=None):
         self.error_msg = """\nDW driver accepts 0 or 4 arguments.\nExample: python driver.py -m DoubleWell -o omega_b (cm^-1) V0 (cm^-1) mass(a.u) delta(angs) \n
         python driver.py -m DoubleWell -o 500,2085,1837,0.00 \n"""
-        super(DoubleWell_driver, self).__init__(args)
+        super(DoubleWell_driver, self).__init__(args, error_msg=self.error_msg)
 
     def check_arguments(self):
         """Function that checks the arguments required to run the driver"""
@@ -46,8 +44,7 @@ def check_arguments(self):
             self.delta = 00
         else:
             try:
-                arglist = self.args.split(",")
-                param = list(map(float, arglist))
+                param = list(map(float, self.args))
                 assert len(param) == 4
                 w_b = param[0] * invcm2au
                 v0 = param[1] * invcm2au

drivers/py/pes/doublewell_with_bath.py

@@ -16,7 +16,6 @@
 __DRIVER_NAME__ = "DW_bath"
 __DRIVER_CLASS__ = "DoubleWell_with_explicit_bath_driver"
 
-# np.set_printoptions(precision=14, suppress=True,threshold='nan',linewidth=1000)
 
 invcm2au = units.unit_to_internal("frequency", "inversecm", 1.0)
 A2au = units.unit_to_internal("length", "angstrom", 1.0)
@@ -44,16 +43,17 @@ def __init__(self, args=None, verbose=None):
         self.error_msg = """\nDW+explicit_bath driver expects 9 arguments.\n
         Example: python driver.py -m DoubleWell_with_explicit_bath -o omega_b (cm^-1) V0 (cm^-1) mass delta(\AA) eta0  eps1 eps2  deltaQ omega_c(cm^-1)     \n
         python driver.py -m DoubleWell -o 500,2085,1837,0.00,1,0,0,1,500\n"""
-        self.args = args
-        self.check_arguments()
+        super(DoubleWell_with_explicit_bath_driver, self).__init__(
+            args, error_msg=self.error_msg
+        )
+
         self.init = False
 
     def check_arguments(self):
         """Function that checks the arguments required to run the driver"""
 
         try:
-            arglist = self.args.split(",")
-            param = list(map(float, arglist))
+            param = list(map(float, self.args))
             assert len(param) == 9
             w_b = param[0] * invcm2au
             v0 = param[1] * invcm2au
@@ -62,7 +62,7 @@ def check_arguments(self):
 
             self.bath_parameters = {}
             self.bath_parameters["m"] = param[2]
-            # self.bath_parameters["delta"] = param[3] * A2au
+
             self.bath_parameters["eta0"] = param[4]
             self.bath_parameters["eps1"] = param[5]
             self.bath_parameters["eps2"] = param[6]

drivers/py/pes/doublewell_with_friction.py

@@ -15,7 +15,6 @@
 __DRIVER_NAME__ = "DW_friction"
 __DRIVER_CLASS__ = "DoubleWell_with_friction_driver"
 
-# np.set_printoptions(precision=14, suppress=True,threshold='nan',linewidth=1000)
 
 invcm2au = units.unit_to_internal("frequency", "inversecm", 1.0)
 A2au = units.unit_to_internal("length", "angstrom", 1.0)
@@ -42,16 +41,16 @@ def __init__(self, args=None, verbose=None):
         self.error_msg = """\nDW+fric driver expects 8 arguments.\n
         Example: python driver.py -m DoubleWell_with_fric -o omega_b (cm^-1) V0 (cm^-1) mass delta(\AA) eta0  eps1 eps2  deltaQ      \n
         python driver.py -m DoubleWell -o 500,2085,1837,0.00,1,0,0,1\n"""
-        self.args = args
+        self.args = args.split(",")
+        self.verbose = verbose
         self.check_arguments()
 
     def check_arguments(self):
         """Function that checks the arguments required to run the driver"""
 
         self.k = 1837.36223469 * (3800.0 / 219323.0) ** 2
         try:
-            arglist = self.args.split(",")
-            param = list(map(float, arglist))
+            param = list(map(float, self.args))
             assert len(param) == 8
             w_b = param[0] * invcm2au
             v0 = param[1] * invcm2au