Fixed some bugs in the new PRNG (basically direct calls to prng.rng f…

…ail) Also created examples of atomswap and alchemical swaps, so that at least they get tested for running properly
lab-cosmo · Apr 29, 2024 · fe94cd4 · fe94cd4
1 parent 8fdf525
commit fe94cd4
Show file tree

Hide file tree

Showing 13 changed files with 373 additions and 17 deletions.
diff --git a/drivers/f90/LJ.f90 b/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
diff --git a/drivers/f90/driver.f90 b/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

diff --git a/examples/features/monte_carlo/alchemical/init.xyz b/examples/features/monte_carlo/alchemical/init.xyz
@@ -0,0 +1,8 @@
+6
+# CELL(abcABC):  130.18034   170.29206   240.21758    90.00000    90.00000    90.00000  cell{atomic_unit}  Traj: positions{angstrom} Step:          12  Bead:       0 
+H 1.3762859974402626 1.6475377505750781 1.3947309252535616  
+C 0.5610492771673459 0.6910172235578630 0.5785825986233499  
+H 0.0313332395759015 0.1474340996450444 1.3607366615876173  
+D 1.3412729333952838 0.1418736998062487 0.0516893010969587  
+D -0.0352236727959057 1.3786671745683925 -0.0209976311964991  
+H -0.162629  -0.151689  -0.290930
diff --git a/examples/features/monte_carlo/alchemical/input.xml b/examples/features/monte_carlo/alchemical/input.xml
@@ -0,0 +1,42 @@
+<simulation verbosity='medium'>
+  <output prefix='simulation'>
+    <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, pressure_cv{megapascal} ] </properties>
+    <trajectory filename='pos' stride='20' cell_units='angstrom' format="xyz"> positions{angstrom} </trajectory>
+    <checkpoint stride='200'/>
+  </output>
+  <total_steps>1000</total_steps>
+  <prng>
+    <seed>31415</seed>
+  </prng>
+  <ffsocket name='ch5p' mode='unix' pbc='false'>
+      <address>ch5p</address>
+  </ffsocket>
+  <system>
+    <initialize nbeads='4'>
+      <file mode='xyz'> init.xyz </file>
+      <velocities mode='thermal' units='kelvin'> 100 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='ch5p'> </force>
+    </forces>
+    <motion mode='multi'>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 0.5 </timestep>
+        <thermostat mode='langevin'>
+          <tau units='femtosecond'> 10 </tau>
+        </thermostat>
+      </dynamics>
+    </motion>
+    <motion mode='alchemy'>
+      <alchemy>
+        <nxc> 1 </nxc>
+        <names> [ H, D ] </names>
+      </alchemy>
+    </motion>
+    </motion>
+    <ensemble>
+      <temperature units='kelvin'> 400 </temperature>
+    </ensemble>
+  </system>
+</simulation>
diff --git a/examples/features/monte_carlo/alchemical/run.sh b/examples/features/monte_carlo/alchemical/run.sh
@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m ch4hcbe -u -a ch5p"
+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}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"
diff --git a/examples/features/monte_carlo/atomswap/input.xml b/examples/features/monte_carlo/atomswap/input.xml
@@ -0,0 +1,42 @@
+<simulation verbosity='high'>
+  <output prefix='simulation'>
+    <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, pressure_cv{megapascal} ] </properties>
+    <trajectory filename='pos' stride='20' cell_units='angstrom' format="ase"> positions{angstrom} </trajectory>
+    <checkpoint stride='200'/>
+  </output>
+  <total_steps>1000</total_steps>
+  <prng>
+    <seed>31415</seed>
+  </prng>
+  <ffsocket name='ljmix' mode='unix' pbc='false'>
+      <address>ljmix</address>
+  </ffsocket>
+  <system>
+    <initialize nbeads='1'>
+      <file mode='xyz'> ljmix.xyz </file>
+      <velocities mode='thermal' units='kelvin'> 100 </velocities>
+    </initialize>
+    <forces>
+      <force forcefield='ljmix'> </force>
+    </forces>
+    <motion mode='multi'>
+    <motion mode='dynamics'>
+      <dynamics mode='nvt'>
+        <timestep units='femtosecond'> 5 </timestep>
+        <thermostat mode='svr'>
+          <tau units='femtosecond'> 10 </tau>
+        </thermostat>
+      </dynamics>
+    </motion>
+    <motion mode='atomswap'>
+      <atomswap>
+        <nxc> 1 </nxc>
+        <names> [ Ne, Ar ] </names>
+      </atomswap>
+    </motion>
+    </motion>
+    <ensemble>
+      <temperature units='kelvin'> 300 </temperature>
+    </ensemble>
+  </system>
+</simulation>
diff --git a/examples/features/monte_carlo/atomswap/ljmix.xyz b/examples/features/monte_carlo/atomswap/ljmix.xyz
@@ -0,0 +1 @@
+../../../init_files/ljmix.xyz
diff --git a/examples/features/monte_carlo/atomswap/run.sh b/examples/features/monte_carlo/atomswap/run.sh
@@ -0,0 +1,16 @@
+ipi=i-pi
+driver="i-pi-driver -m ljmix -u -a ljmix -o 32,4,1e-3,12"
+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}
+
+${driver} > /dev/null & 
+echo "# Driver is running"
+
+wait
+
+echo "# Simulation complete"
diff --git a/examples/init_files/ljmix.xyz b/examples/init_files/ljmix.xyz
@@ -0,0 +1,110 @@
+108
+#  CELL(abcABC):  17.8469984966  17.8469984966  17.8469984966  90.00000  90.00000  90.00000  Traj:  positions{angstrom}  Step:  0  Bead:  0  cell{angstrom}
+Ne    0.147000137166   0.132000079971   0.114000117173
+Ne    -0.125000123837   2.94299788633   3.18800164278
+Ne    3.04200165055   2.93700230854   0.107000161039
+Ne    2.97599737714   -0.147999752916   2.96000035008
+Ne    -0.228999850096   -0.0300000085539   5.92297467609
+Ne    0.0939998645211   2.99900071046   8.91499553915
+Ne    2.95399948052   2.89500151338   5.91000983444
+Ne    2.8920010786   -0.155000238227   8.8700154763
+Ne    -0.0480000242697   0.0869999083873   11.900983782
+Ne    0.0169999766244   2.94699846603   14.9550242141
+Ne    2.97899781192   2.98899926119   11.875001181
+Ne    2.97599737714   0.114000117173   14.7939955891
+Ne    -0.0529997434675   5.85201201223   -0.165000099964
+Ne    0.00600002287786   8.85800315363   3.04000136069
+Ne    3.04100150562   8.95097958943   -0.0960001543749
+Ne    3.00700186988   5.87900004994   2.87599875977
+Ne    0.0569997939979   5.91000983444   6.03600692814
+Ne    -0.0969997701246   8.8669991662   8.80397416049
+Ne    2.83199767476   8.76597923681   5.71897686163
+Ne    3.14800113748   5.95599533399   8.96399734879
+Ne    0.0309999947276   5.97398735913   11.9319935665
+Ne    0.12399997891   8.97198792467   14.7229800075
+Ne    2.99800056553   8.86101946373   11.9519964651
+Ne    2.97100194428   5.94297757463   14.9140129803
+Ne    0.00899998139907   11.7700124225   0.0890001982411
+Ne    0.0150000042769   14.9309995688   2.99199969597
+Ne    2.92400042449   14.8889828983   0.0619999894552
+Ne    2.83199767476   11.9180232882   3.16199787469
+Ne    -0.0969997701246   11.8850026303   6.10501163633
+Ne    0.169000150494   14.7310235011   8.88901293814
+Ne    2.8769989047   14.7500209629   5.92101672041
+Ne    3.07800157614   12.0379877617   8.93600387439
+Ar   0.0190000018896   12.010999724   11.8850026303
+Ar   0.139000036105   14.8480245822   14.8409865253
+Ar   2.90399752595   14.8350068229   12.0199957365
+Ar   3.08300230078   12.0130105974   14.9420064547
+Ar   6.04198663062   0.050999982791   -0.122000218234
+Ar   5.91101527114   2.9569999153   2.81099992306
+Ar   8.94902163375   3.11100106696   0.0129999790117
+Ar   8.9909853865   -0.0100000204905   2.79399745931
+Ar   6.06098409246   0.0420000013919   6.018014903
+Ar   6.15602431937   3.00200114524   8.7079814146
+Ar   8.92198067832   3.0690002718   5.90900439774
+Ar   8.98802199413   0.132000079971   9.00601401927
+Ar   5.91598953692   0.156999998903   11.9360153133
+Ar   6.01097684611   2.8409989791   14.9089857968
+Ar   8.97399879806   2.98999940612   11.9060109655
+Ar   8.94600532365   -0.011999992838   14.7619803679
+Ar   5.87100947407   5.8559808413   -0.118999783453
+Ar   6.16401489524   8.84297452087   3.01699802737
+Ar   9.19302524528   8.76201040774   0.075999901723
+Ar   8.8339785083   5.96499134656   2.83499810954
+Ar   5.9870051185   6.00298627024   5.94101961895
+Ar   5.97700366923   8.97399879806   8.83900569179
+Ar   9.02501148111   8.84599083096   6.02198373208
+Ar   8.84800170436   5.69399969732   9.14301799894
+Ar   5.92398011279   6.03198518135   11.8649997317
+Ar   6.0190203397   8.82699336913   14.8660166074
+Ar   8.834983945   9.07602416415   11.9409895791
+Ar   8.89102381154   6.03902323824   14.8950155185
+Ar   5.8559808413   11.9979819646   0.104999871185
+Ar   5.82798736689   14.7829887031   2.90299738102
+Ar   8.94198357686   14.8989843475   -0.0420000013919
+Ar   8.93002417191   11.9660196611   3.01499773752
+Ar   6.07299641512   11.8649997317   5.91202070784
+Ar   6.03097974465   14.8909937717   8.97902598156
+Ar   8.97198792467   14.8580260315   5.88302179673
+Ar   8.95701220962   11.8349953839   8.96701365889
+Ar   5.96599678326   11.8819863202   12.0440203819
+Ar   5.94699932142   14.7049879824   14.7739926905
+Ar   8.87599517877   14.852998848   11.9509910284
+Ar   8.9919908232   11.9170178515   14.9289886954
+Ar   11.9460167626   0.0880000533142   -0.0279999832887
+Ar   11.8290156815   3.08100201092   3.11700193652
+Ar   15.0040260237   3.0330003462   -0.00799999522534
+Ar   14.8280216837   0.11999992838   3.0289997665
+Ar   11.9060109655   0.0169999766244   5.99499569437
+Ar   11.8729903076   2.91599926507   9.01601546854
+Ar   14.9069749234   3.06199925731   5.96001708079
+Ar   14.9479861572   -0.122000218234   8.8639828561
+Ar   11.8690214785   -0.118000167703   11.9299826931
+Ar   11.9319935665   2.86300216749   14.9789959417
+Ar   14.8350068229   2.97599737714   12.0179848631
+Ar   14.8519934113   -0.217999843432   14.7829887031
+Ar   11.8919877694   5.99002142859   0.0579999389248
+Ar   11.8159979221   8.68400968698   2.93600216361
+Ar   14.9579876065   8.87599517877   -0.0550000333213
+Ar   14.8399810886   6.06098409246   2.91599926507
+Ar   11.9109852313   5.97102396676   5.72601491853
+Ar   11.9209866806   8.94701076035   8.91499553915
+Ar   15.0320194981   8.97198792467   5.83100367699
+Ar   14.6990082799   5.88900149921   8.93097669089
+Ar   11.764985239   5.84402143636   11.8460022699
+Ar   11.9830062496   8.75502526857   14.9330104422
+Ar   14.8399810886   8.92600242512   11.9470221993
+Ar   14.9759796316   5.71199172246   14.7709763805
+Ar   11.823988498   11.8709794342   -0.0070000090516
+Ar   11.8939986428   14.8399810886   2.77099941776
+Ar   15.0390046373   14.9279832587   0.0209999742371
+Ar   14.8760180566   11.8509765357   2.90999839551
+Ar   11.7849881376   11.9340044399   6.07098554172
+Ar   11.93702075   14.6810162548   8.92901873521
+Ar   14.9350213156   14.9020006576   5.87302034746
+Ar   14.8819977591   11.9440058892   9.00701945597
+Ar   11.9409895791   11.9560182119   11.8770120544
+Ar   11.7910207578   14.8999897842   14.9550242141
+Ar   14.838023133   14.9469807205   11.8309736371
+Ar   14.8909937717   12.0050200215   15.0500115233