Merge pull request #53 from ltimmerman3/master

Adding examples for run_mode_b for simulations of type singlepoint, r…

examples/FileIO/md/ab-initio/run.py

@@ -0,0 +1,32 @@
+from sparc.calculator import SPARC
+from ase import Atoms
+import numpy as np
+
+Ag_cluster = Atoms('Ag5', positions = [(0.0, 2.6579, 0.9366), (0.0, -1.3587, -1.4045), (0.0, 0.0, 0.9358), (0.0, -2.6579, 0.9366),
+                                      (0.0, 1.3587, -1.4045)],
+                  pbc = (0,0,0))
+Ag_cluster.set_cell([20., 24., 24.])
+Ag_cluster.center()
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+    "MD_FLAG": 1,
+    "MD_METHOD": "NVK_G",
+    "ION_TEMP": 10,
+    "MD_NSTEP": 10,
+    "MD_TIMESTEP": 2,
+}
+
+Ag_cluster.calc = SPARC(label = "Ag_cluster", **calc_params)
+energy = Ag_cluster.get_potential_energy()
+print("Energy: ", energy)
+forces = Ag_cluster.get_forces()
+print("Max Force: ", np.max(forces))

examples/FileIO/md/mlff/run.py

@@ -0,0 +1,34 @@
+from sparc.calculator import SPARC
+from ase import Atoms
+import numpy as np
+
+Ag_cluster = Atoms('Ag5', positions = [(0.0, 2.6579, 0.9366), (0.0, -1.3587, -1.4045), (0.0, 0.0, 0.9358), (0.0, -2.6579, 0.9366),
+                                      (0.0, 1.3587, -1.4045)],
+                  pbc = (0,0,0))
+Ag_cluster.set_cell([20., 24., 24.])
+Ag_cluster.center()
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+    "MD_FLAG": 1,
+    "MD_METHOD": "NVK_G",
+    "ION_TEMP": 10,
+    "MD_NSTEP": 10,
+    "MD_TIMESTEP": 2,
+    "MLFF_FLAG": 1,
+    "MLFF_INITIAL_STEPS_TRAIN": 3
+}
+
+Ag_cluster.calc = SPARC(label = "Ag_cluster", **calc_params)
+energy = Ag_cluster.get_potential_energy()
+print("Energy: ", energy)
+forces = Ag_cluster.get_forces()
+print("Max Force: ", np.max(forces))

examples/FileIO/relax/relax_cell/ab-initio/run.py

@@ -0,0 +1,33 @@
+# Using SPARC
+from sparc.calculator import SPARC
+from ase.build import bulk
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [4,4,4],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "CALC_STRESS": 1,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+    "RELAX_FLAG": 2,
+}
+
+
+atoms = bulk('Rh', crystalstructure='fcc', a = 3.81)
+atoms.calc = SPARC(**calc_params)
+# Trigger the calculation by calling a property that will be available after the calculation
+stress = atoms.get_stress()
+# Get the final energy and volume from a single point calculation
+del calc_params['RELAX_FLAG']
+if type(atoms) == list:
+    atoms = atoms[-1]
+atoms.calc = SPARC(**calc_params)
+energy = atoms.get_potential_energy()
+volume = atoms.get_volume()
+
+print('***********************************************************************************************')
+print('v0 = {0} A^3\nE0 = {1} eV'.format(volume, energy))

examples/FileIO/relax/relax_coords/ab-initio/run.py

@@ -0,0 +1,29 @@
+# Using SPARC
+from sparc.calculator import SPARC
+from ase.build import molecule
+import numpy as np
+
+water = molecule('H2O', vacuum=7)
+water.pbc = [False,False,False]
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "ELEC_TEMP_TYPE": "fermi-dirac",
+    "ELEC_TEMP": 116,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+    "RELAX_FLAG": 1,
+}
+
+water.calc = SPARC(**calc_params)
+energy = water.get_potential_energy()
+forces = water.get_forces()
+
+print('Energy:', energy)
+print('Max Force:', np.max(forces))

examples/FileIO/relax/relax_coords/mlff/run.py

@@ -0,0 +1,31 @@
+# Using SPARC
+from sparc.calculator import SPARC
+from ase.build import molecule
+import numpy as np
+
+water = molecule('H2O', vacuum=7)
+water.pbc = [False,False,False]
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "ELEC_TEMP_TYPE": "fermi-dirac",
+    "ELEC_TEMP": 116,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+    "RELAX_FLAG": 1,
+    "MLFF_FLAG": 1,
+    "MLFF_INITIAL_STEPS_TRAIN": 3,
+}
+
+water.calc = SPARC(**calc_params)
+energy = water.get_potential_energy()
+forces = water.get_forces()
+
+print('Energy:', energy)
+print('Max Force:', np.max(forces))

examples/env.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+# set path to any shared libraries such as mkl, blas, scalapack etc.
+ml intel-oneapi-compilers intel-oneapi-mkl intel-oneapi-mpi
+source ~/p-amedford6-0/venvs/joss-sparc-x-api/bin/activate
+
+cd ~
+export HOME=$(pwd)
+cd -
+export SCRATCH=$HOME/scratch
+export DATA=$HOME/p-amedford6-0
+
+export PYTHONPATH=$DATA/active/JOSS-SPARC-X-API/SPARC-X-API:$PYTHONPATH
+
+# Check if Slurm is installed and running
+if command -v sinfo &> /dev/null; then
+  export ASE_SPARC_COMMAND="srun ${DATA}/active/JOSS-SPARC-X-API/dev_SPARC/lib/sparc --export=ALL -K"
+else
+  nprocs=$(nproc)
+  export ASE_SPARC_COMMAND="mpirun -np $nprocs ${DATA}/SPARC/lib/sparc --export=ALL -K" 
+fi
+# Set pseudopotential path is you did not download with python -m sparc.download_data
+export PLUMED_KERNEL=/storage/home/hcoda1/9/ltimmerman3/p-amedford6-0/active/JOSS-SPARC-X-API/plumed-2.9.2/src/lib/libplumedKernel.so

examples/socket/md/ab-initio/ase/run.py

@@ -0,0 +1,35 @@
+from sparc.calculator import SPARC
+from ase import Atoms
+from ase.md.nvtberendsen import NVTBerendsen
+from ase.md.verlet import VelocityVerlet
+from ase.units import fs
+from ase.constraints import FixedPlane
+
+Ag_cluster = Atoms('Ag5', positions = [(0.0, 2.6579, 0.9366), (0.0, -1.3587, -1.4045), (0.0, 0.0, 0.9358), (0.0, -2.6579, 0.9366),
+                                      (0.0, 1.3587, -1.4045)],
+                  pbc = (0,0,0))
+Ag_cluster.set_cell([20., 24., 24.])
+Ag_cluster.center()
+cons = [FixedPlane(i, [1, 0, 0]) for i in range(5)]
+Ag_cluster.set_constraint(cons)
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+}
+
+taut = 50 * fs
+timestep = 2 * fs
+with SPARC(use_socket=True, **calc_params) as calc:
+    Ag_cluster.calc = calc
+    #dyn = NVTBerendsen(Ag_cluster, temperature_K=10, taut=taut, timestep=timestep, trajectory='Ag-cluster-md.traj')
+    dyn = VelocityVerlet(Ag_cluster, dt=2.0 * fs,
+                     trajectory='Ag-cluster-verlet.traj', logfile='verlet-md.log')
+    dyn.run(10)

examples/socket/md/ab-initio/plumed/run.py

@@ -0,0 +1,84 @@
+# Using SPARC
+from sparc.calculator import SPARC
+from ase import units
+from ase import Atoms
+from ase.calculators.plumed import Plumed
+from ase.md.nvtberendsen import NVTBerendsen
+
+calc_params = {
+    "EXCHANGE_CORRELATION": "GGA_PBE",
+    "KPOINT_GRID": [1,1,1],
+    "MESH_SPACING": 0.35,
+    "TOL_SCF": 0.0001,
+    "MAXIT_SCF": 100,
+    "PRINT_RESTART_FQ": 10,
+    "PRINT_ATOMS": 1,
+    "PRINT_FORCES": 1,
+    "SPIN_TYP": 0,
+}
+
+timestep = 2 * units.fs  # the units module contains the conversion factor to go from units.<unit> to ASE units via multiplication
+ps = 1000 * units.fs
+
+# Units conversion: everything will be in ASE units, so a time step is in ASE unit and needs to be interpreted as some multiple of sensical units
+# Likewise, energy is in eV which is ~98.6 kJ/mol, so a single unit of energy is 98.6 kJ/mol which are the internal units of PLUMED
+setup = [
+    f"UNITS LENGTH=A TIME={1/ps} ENERGY={units.mol/units.kJ}",  # Set units to match desired properties
+
+    # Calculate the center of mass of atoms 1-5
+    "com: COM ATOMS=1-5",
+
+    # Define the coordination number (C)
+    "c: COORDINATION GROUPA=1-5 SWITCH={RATIONAL R_0=3.0 NN=8 MM=16} NOPBC",
+
+    # Define the radius of gyration (R)
+    "r: GYRATION TYPE=RADIUS ATOMS=1-5 NOPBC",
+
+    # Compute CV1 and CV2 as orthogonal linear combinations of C and R
+    "cv1: COMBINE ARG=c,r COEFFICIENTS=0.99715,-0.07534 PERIODIC=NO",
+    "cv2: COMBINE ARG=c,r COEFFICIENTS=0.07534,0.99715 PERIODIC=NO",
+
+    # Apply lower wall on CV1 at 5.0 with harmonic constant 10 eV
+    f"LOWER_WALLS ARG=cv1 AT=5.0 KAPPA=10.0",
+
+    # Apply lower wall on CV2 at 3.0 with harmonic constant 50 eV
+    f"UPPER_WALLS ARG=cv2 AT=3.0 KAPPA=50.0",
+
+    # Perform well-tempered metadynamics on CV1 and CV2
+    f"METAD ARG=cv1,cv2 HEIGHT=0.3 PACE=500 SIGMA=0.3,0.03 GRID_MIN=0.0,0.0 GRID_MAX=10.0,5.0 GRID_BIN=500,500 BIASFACTOR=100 FILE=HILLS",
+
+    # Print out the collective variables for monitoring
+    "PRINT ARG=cv1,cv2,c,r FILE=COLVAR STRIDE=1"
+]
+
+Ag_cluster = Atoms('Ag5', positions = [(0.0, 2.6579, 0.9366), (0.0, -1.3587, -1.4045), (0.0, 0.0, 0.9358), (0.0, -2.6579, 0.9366),
+                                      (0.0, 1.3587, -1.4045)],
+                  pbc = (0,0,0))
+Ag_cluster.set_cell([20., 24., 24.])
+Ag_cluster.center()
+# cons = [FixedPlane(i, [1, 0, 0]) for i in range(5)]
+# Ag_cluster.set_constraint(cons)
+with SPARC(use_socket=True, **calc_params) as calc:
+    Ag_cluster.calc = Plumed(calc=calc,
+                    input=setup,
+                    timestep=timestep,
+                    atoms=Ag_cluster,
+                    kT=0.00861733) # 10 K in eV thermal energy units
+    dyn = NVTBerendsen(Ag_cluster, timestep, temperature_K=0.00861733/units.kB, taut=50*units.fs,
+                fixcm=False, trajectory='Ag-cluster-metadynamics.traj')
+    dyn.run(10)
+
+
+"""
+    # Restrain atoms 1-5 within 2.0 Å of the center of mass using upper walls
+    "d1: DISTANCE ATOMS=1,com",
+    "UPPER_WALLS ARG=d1 AT=2.0 KAPPA=100.0",
+    "d2: DISTANCE ATOMS=2,com",
+    "UPPER_WALLS ARG=d2 AT=2.0 KAPPA=100.0",
+    "d3: DISTANCE ATOMS=3,com",
+    "UPPER_WALLS ARG=d3 AT=2.0 KAPPA=100.0",
+    "d4: DISTANCE ATOMS=4,com",
+    "UPPER_WALLS ARG=d4 AT=2.0 KAPPA=100.0",
+    "d5: DISTANCE ATOMS=5,com",
+    "UPPER_WALLS ARG=d5 AT=2.0 KAPPA=100.0",
+"""