Add spectral abf example to lammps (#332)

examples/lammps/spectral_abf/adp.input

@@ -0,0 +1,244 @@
+LAMMPS data file for ACE
+
+22 atoms
+21 bonds
+36 angles
+66 dihedrals
+
+7 atom types
+8 bond types
+16 angle types
+19 dihedral types
+
+-100 100  xlo xhi
+-100 100  ylo yhi
+-100 100  zlo zhi
+
+Masses
+
+1 1.008
+2 12.01
+3 12.01
+4 16.0
+5 14.01
+6 1.008
+7 1.008
+
+Atoms
+
+ 1 1 1 0.1123 2.0177833000000001 1.0353295 0.054216800000000002
+ 2 1 2 -0.36620000000000003 1.9998066999999999 2.1250665999999998 0.0538082
+ 3 1 1 0.1123 1.4254165999999999 2.4843904000000001 0.90724000000000005
+ 4 1 1 0.1123 1.5326313 2.4823206 -0.86386980000000002
+ 5 1 3 0.5972000021951126 3.4222907 2.6597886000000002 0.13890169999999999
+ 6 1 4 -0.56790000164633436 4.3607215999999998 1.8984493 0.36315979999999998
+ 7 1 5 -0.41570000000000001 3.5879191000000001 3.9603571 -0.1138393
+ 8 1 6 0.27190000000000003 2.7859197999999998 4.5436604999999997 -0.3046603
+ 9 1 2 0.033699999999999994 4.8919525999999998 4.5905944999999999 -0.067514500000000005
+10 1 7 0.082299999999999998 5.4432974999999999 4.2337999000000002 0.80264219999999997
+11 1 2 -0.1825 5.6616692000000004 4.2349395999999997 -1.3354246999999999
+12 1 1 0.060299999999999992 5.1089564000000003 4.5854832999999999 -2.2068146999999998
+13 1 1 0.060299999999999992 6.6419078000000003 4.7111339000000001 -1.3103305999999999
+14 1 1 0.060299999999999992 5.7851258000000003 3.1535202 -1.3931851
+15 1 3 0.59730000054877808 4.7222410999999997 6.1004128 0.021521100000000001
+16 1 4 -0.56790000164633436 3.6035081999999998 6.6057630999999999 -0.039872400000000002
+17 1 5 -0.41570000000000001 5.8431179000000002 6.8242649999999996 0.073516899999999996
+18 1 6 0.27190000000000003 6.7415383999999996 6.3635967999999998 0.087583599999999998
+19 1 2 -0.14899999999999999 5.8209701000000003 8.2733497000000007 0.056992800000000003
+20 1 7 0.097600000000000006 4.7879085999999997 8.6215931000000001 0.0518138
+21 1 7 0.097600000000000006 6.3303031000000001 8.6556759000000003 0.94173300000000004
+22 1 7 0.097600000000000006 6.3290771000000001 8.6329496999999993 -0.83830830000000001
+
+Bonds
+
+ 1 3 2   3
+ 2 3 2   4
+ 3 3 1   2
+ 4 3 11 12
+ 5 3 11 13
+ 6 3 11 14
+ 7 5 9  10
+ 8 7 7   8
+ 9 5 19 20
+10 5 19 21
+11 5 19 22
+12 7 17 18
+13 1 5   6
+14 2 5   7
+15 4 2   5
+16 1 15 16
+17 2 15 17
+18 6 9  11
+19 4 9  15
+20 8 7   9
+21 8 17 19
+
+Angles
+
+1 2 5 7 8
+2 4 4 2 5
+3 5 3 2 4
+4 4 3 2 5
+5 5 1 2 3
+6 5 1 2 4
+7 4 1 2 5
+8 2 15 17 18
+9 5 13 11 14
+10 5 12 11 13
+11 5 12 11 14
+12 9 10 9 11
+13 10 10 9 15
+14 11 9 11 12
+15 11 9 11 13
+16 11 9 11 14
+17 12 8 7 9
+18 13 7 9 10
+19 16 21 19 22
+20 16 20 19 21
+21 16 20 19 22
+22 12 18 17 19
+23 13 17 19 20
+24 13 17 19 21
+25 13 17 19 22
+26 1 6 5 7
+27 3 5 7 9
+28 6 2 5 6
+29 7 2 5 7
+30 1 16 15 17
+31 3 15 17 19
+32 8 11 9 15
+33 6 9 15 16
+34 7 9 15 17
+35 14 7 9 11
+36 15 7 9 15
+
+Dihedrals
+
+1 1 6 5 7 8
+2 2 6 5 7 8
+3 3 5 7 9 10
+4 10 4 2 5 6
+5 3 4 2 5 6
+6 11 4 2 5 6
+7 3 4 2 5 7
+8 10 3 2 5 6
+9 3 3 2 5 6
+10 11 3 2 5 6
+11 3 3 2 5 7
+12 2 2 5 7 8
+13 10 1 2 5 6
+14 3 1 2 5 6
+15 11 1 2 5 6
+16 3 1 2 5 7
+17 1 16 15 17 18
+18 2 16 15 17 18
+19 3 15 17 19 20
+20 3 15 17 19 21
+21 3 15 17 19 22
+22 12 14 11 9 15
+23 12 13 11 9 15
+24 12 12 11 9 15
+25 12 10 9 11 12
+26 12 10 9 11 13
+27 12 10 9 11 14
+28 10 10 9 15 16
+29 11 10 9 15 16
+30 3 10 9 15 17
+31 2 9 15 17 18
+32 3 8 7 9 10
+33 3 8 7 9 11
+34 3 8 7 9 15
+35 12 7 9 11 12
+36 12 7 9 11 13
+37 12 7 9 11 14
+38 3 18 17 19 20
+39 3 18 17 19 21
+40 3 18 17 19 22
+41 19 5 9 7 8
+42 19 15 19 17 18
+43 2 6 5 7 9
+44 1 5 7 9 11
+45 4 5 7 9 11
+46 5 5 7 9 11
+47 6 5 7 9 11
+48 7 5 7 9 15
+49 8 5 7 9 15
+50 9 5 7 9 15
+51 6 5 7 9 15
+52 2 2 5 7 9
+53 2 16 15 17 19
+54 3 11 9 15 16
+55 13 11 9 15 17
+56 14 11 9 15 17
+57 5 11 9 15 17
+58 6 11 9 15 17
+59 2 9 15 17 19
+60 3 7 9 15 16
+61 15 7 9 15 17
+62 16 7 9 15 17
+63 17 7 9 15 17
+64 6 7 9 15 17
+65 18 2 7 5 6
+66 18 9 17 15 16
+
+Pair Coeffs
+
+1 0.01570000002623629  2.6495327872602221
+2 0.10939999991572773  3.3996695084507409
+3 0.086000000128358844 3.3996695079448309
+4 0.20999999984182244  2.9599219016446874
+5 0.16999999991766696  3.2499985240310356
+6 0.015700000004219245 1.0690784617205229
+7 0.015700000098461422 2.4713530426421655
+
+Bond Coeffs
+
+1 570.0 1.229
+2 490.0 1.335
+3 340.0 1.090
+4 317.0 1.522
+5 340.0 1.090
+6 310.0 1.526
+7 434.0 1.010
+8 337.0 1.449
+
+Angle Coeffs
+
+1 80.0 122.90005267195104
+2 50.0 120.00005142908158
+3 50.0 121.90005224337536
+4 50.0 109.50004692903693
+5 35.0 109.50004692903693
+6 80.0 120.40005160051184
+7 70.0 116.60004997192425
+8 63.0 111.10004761475803
+9 50.0 109.50004692903693
+10 50.0 109.50004692903693
+11 50.0 109.50004692903693
+12 50.0 118.04005047448166
+13 50.0 109.50004692903693
+14 80.0 109.70004701475206
+15 63.0 110.10004718618234
+16 35.0 109.50004692903693
+
+Dihedral Coeffs
+
+ 1 2.0          1 1  
+ 2 2.5         -1 2  
+ 3 0.0          1 2  
+ 4 2.0          1 2  
+ 5 0.4          1 3  
+ 6 0.0          1 4  
+ 7 0.0          1 1  
+ 8 0.272        1 2  
+ 9 0.43         1 3  
+10 0.8          1 1  
+11 0.08        -1 3  
+12 0.155555556  1 3  
+13 0.20         1 1  
+14 0.20         1 2  
+15 0.45        -1 1  
+16 1.58        -1 2  
+17 0.55        -1 3  
+18 10.5        -1 2  
+19 1.10        -1 2  

examples/lammps/spectral_abf/adp.lmp

@@ -0,0 +1,24 @@
+units			real
+
+neigh_modify    once yes  one  22 page 2200
+
+atom_style	full
+bond_style      harmonic
+angle_style     harmonic
+dihedral_style  harmonic
+pair_style      lj/cut/coul/cut 10.0
+pair_modify     mix arithmetic
+
+kspace_style    none
+read_data       adp.input
+
+
+velocity all create 300 3
+
+timestep        1.0
+
+fix             1 all nve
+fix             2 all langevin 300 300 1000 63683
+fix             3 all momentum 100 linear 1 1 1
+fix             4 all shake 0.0001 10 100 b 3 5 7
+special_bonds   amber

examples/lammps/spectral_abf/run.py

@@ -0,0 +1,107 @@
+#!/usr/bin/env python3
+
+"""
+Example SpectralABF simulation with pysages and lammps.
+
+For a list of possible options for running the script pass `-h` as argument from the
+command line, or call `get_args(["-h"])` if the module was loaded interactively.
+"""
+
+# %%
+import argparse
+import sys
+
+import numpy
+from lammps import lammps
+
+import pysages
+from pysages.colvars import DihedralAngle
+from pysages.methods import HistogramLogger, SpectralABF
+
+
+# %%
+def generate_context(args="", script="adp.lmp", store_freq=1):
+    """
+    Returns a lammps simulation defined by the contents of `script` using `args` as
+    initialization arguments.
+    """
+    context = lammps(cmdargs=args.split())
+    context.file(script)
+    # Allow for the retrieval of the wrapped positions
+    context.command(f"dump 4a all custom {store_freq} dump.myforce id type x y z")
+    return context
+
+
+def get_args(argv):
+    """Process the command-line arguments to this script."""
+
+    available_args = [
+        ("time-steps", "t", int, 2e6, "Number of simulation steps"),
+        ("kokkos", "k", bool, True, "Whether to use Kokkos acceleration"),
+        ("log-steps", "l", int, 2e3, "Number of simulation steps for logging"),
+    ]
+    parser = argparse.ArgumentParser(description="Example script to run pysages with lammps")
+
+    for name, short, T, val, doc in available_args:
+        if T is bool:
+            action = "store_" + str(val).lower()
+            parser.add_argument("--" + name, "-" + short, action=action, help=doc)
+        else:
+            convert = (lambda x: int(float(x))) if T is int else T
+            parser.add_argument("--" + name, "-" + short, type=convert, default=T(val), help=doc)
+
+    return parser.parse_args(argv)
+
+
+def main(argv):
+    """Example SpectralABF simulation with pysages and lammps."""
+    args = get_args(argv)
+
+    context_args = {"store_freq": args.log_steps}
+    if args.kokkos:
+        # Passed to the lammps constructor as `cmdargs` when running the script
+        # with the --kokkos (or -k) option
+        context_args["args"] = "-k on g 1 -sf kk -pk kokkos newton on neigh half"
+    #        context_args["args"] = "-k on -sf kk -pk kokkos newton on neigh half"
+
+    # Setting the collective variable, method, and running the simulation
+    cvs = [DihedralAngle([4, 6, 8, 14]), DihedralAngle([6, 8, 14, 16])]
+    grid = pysages.Grid(
+        lower=(-numpy.pi, -numpy.pi),
+        upper=(numpy.pi, numpy.pi),
+        shape=(32, 32),
+        periodic=True,
+    )
+    method = SpectralABF(cvs, grid)
+    callback = HistogramLogger(args.log_steps)
+    raw_result = pysages.run(
+        method,
+        generate_context,
+        args.time_steps,
+        callback=callback,
+        context_args=context_args,
+    )
+    # Post-run analysis
+    result = pysages.analyze(raw_result)
+    mesh = result["mesh"]
+    fes_fn = result["fes_fn"]
+    A = fes_fn(mesh)
+    hist = result["histogram"]
+    A = A.reshape(32, 32)
+    numpy.savetxt("FES.dat", numpy.hstack([mesh, A.reshape(-1, 1)]))
+    numpy.savetxt("hist.dat", numpy.hstack([mesh, hist.reshape(-1, 1)]))
+    bins = 50
+    histo, xedges, yedges = numpy.histogram2d(
+        callback.data[:, 0],
+        callback.data[:, 1],
+        bins=bins,
+        range=[[-numpy.pi, numpy.pi], [-numpy.pi, numpy.pi]],
+    )
+    xedges = (xedges[1:] + xedges[:-1]) / 2
+    yedges = (yedges[1:] + yedges[:-1]) / 2
+    mesh = numpy.dstack(numpy.meshgrid(xedges, yedges)).reshape(-1, 2)
+    numpy.savetxt("hist-from-logger.dat", numpy.hstack([mesh, histo.reshape(-1, 1)]))
+
+
+if __name__ == "__main__":
+    main(sys.argv[1:])