added 1000 ns run files and corrected NPT eqbm time

parameters/mdp/molecular-dynamics/npt_eqbm.mdp

@@ -5,7 +5,7 @@
 ;;; Run control ;;;
 integrator                = md           ; "steep" for energy minimization and "md" for a leap frog algorithm for integrating Newton's eq of motion
 dt                        = 0.001        ; 1 femtosecond
-nsteps                    = 500000      ; 1 nanosecond
+nsteps                    = 5000000      ; 5 nanoseconds
 comm-mode                 = Linear       ; "Linear" (remove center for mass translation) "Angular" ( >> >> and rotation around the center of mass) "None" (no restric) 
 nstcomm                   = 100          ; [steps] frequency for center of mass motion removal
 comm-grps                 = Aqueous Crystal   ; groups for which the center of mass motion is removed

parameters/mdp/molecular-dynamics/npt_prod_1000ns.mdp

@@ -0,0 +1,61 @@
+; Created by Alec Glisman (GitHub: @alec-glisman) on December 31st, 2021
+;;; GROMACS Documentation: https://manual.gromacs.org/documentation/current/user-guide/mdp-options.html ;;;
+;;; Defines ;;;
+; define                   = -DPOSRES_CRYSTAL ; defines a macro that can be used in the topology file
+;;; Run control ;;;
+integrator                = md           ; "steep" for ener minimization and "md" for a leap frog algorithm for integrating Newton's eq of motion
+dt                        = 0.002        ; 2 femtoseconds
+nsteps                    = 500000000    ; 1000 nanoseconds
+comm-mode                 = Linear       ; "Linear" (remove center for mass translation) "Angular" ( >> >> and rotation around the center of mass) "None" (no restric) 
+nstcomm                   = 100          ; [steps] frequency for center of mass motion removal
+comm-grps                 = Aqueous Crystal   ; groups for which the center of mass motion is removed
+;;; Output control ;;;
+nstxout                   = 0           ; [steps] number of steps that elapse between writing coordinates to the output trajectory file (trr)
+nstvout                   = 0           ; [steps] number of steps that elapse between writing velocities to the output trajectory file (trr)
+nstfout                   = 0           ; [steps] number of steps that elapse between writing forces to the output trajectory file (trr)
+nstcalcenergy             = 100         ; [steps] number of steps that elapse between calculating the energies, 0 is never
+nstenergy                 = 1000        ; [steps] number of steps that elapse between writing energies to energy file
+nstlog                    = 1000        ; [steps] number of steps that elapse between writing energies to the log file
+nstxout-compressed        = 1000        ; [steps] number of steps that elapse between writing to .xtc file
+;;; Neighbor searching ;;;
+cutoff-scheme             = Verlet       ; algorithm to generate neigbor list
+nstlist                   = 10           ; [steps] frequency to update the neighbor list
+pbc                       = xyz          ; periodic boundary conditions
+rlist                     = 1.2000       ; [nm] cut-off distance for the short-range neighbor list
+;;; Electrostatics ;;;
+coulombtype               = PME          ; algorithm to generate electrostatics
+rcoulomb                  = 1.2000       ; [nm] The distance for the Coulomb cut-off
+;;; Van der Waals ;;;
+vdw-type                  = cut-off      ; algorithm to generate Van der Waals
+rvdw                      = 1.2000       ; [nm] distance for the LJ or Buckingham cut-off
+DispCorr                  = AllEnerPres  ; corrections to apply for long-ranged energy and/or pressure
+;;; Ewald ;;;
+fourierspacing            = 0.12         ; [nm] for ordinary Ewald, the ratio of the box dimensions and the spacing determines a lower bound for the number of wave vectors to use in each (signed) direction
+pme-order                 = 4            ; interpolation order for PME. 4 equals cubic interpolation. 
+ewald-rtol                = 1.0e-5       ; the relative strength of the Ewald-shifted direct potential at rcoulomb
+ewald-geometry            = 3d           ; "3d" is the ewald sum which is performed in all the three dimensions or "3dc" the reciprocal sum is still performed in 3D but a forec and potential correction applied in the z dimension to produce a pseudo 2D summation
+;;; Temperature coupling ;;;
+Tcoupl                    = nose-hoover            ; Temp coupling using velocity rescaling. Temperature coupling using a Nose-Hoover extended ensemble. The reference temperature and coupling groups are selected as above, but in this case tau-t controls the period of the temperature fluctuations at equilibrium, which is slightly different from a relaxation time. For NVT simulations the conserved energy quantity is written to the energy and log files.
+tc-grps                   = System                 ; groups to couple separately to Temp bath
+tau-t                     = 0.50                   ; [ps] time constant for coupling (one for each group in tc-grps)
+ref-t                     = 300                    ; [K] reference temperature for coupling (one for each group in tc-grps)
+;;; Pressure Coupling ;;;
+Pcoupl                    = Parrinello-Rahman     ; No pressure coupling, This means fixed box size. Extended-ensemble pressure coupling where the box vectors are subject to an equation of motion. The equation of motion for the atoms is coupled to this. No instantaneous scaling takes place.
+Pcoupltype                = semiisotropic         ; specifies the kind of isotropy of the pressure coupling used.
+tau-p                     = 5.0                   ; [ps] The time constant for pressure coupling (one value for all directions).
+ref-p                     = 1.0 1.0               ; [bar] the reference pressure for coupling
+compressibility           = 4.5e-15 4.5e-5        ; [bar^{-1}] for water at 1 atm and 300K the compressibility is 4.5e-5 bar-1
+;;; Velocity generation ;;;
+gen-vel                   = no           ; "no" do not generate velocities "yes" generate velocities in grompp at temp gen-temp which is only meanigful with integrator "md"
+gen-temp                  = 300          ; [K] temperature for Maxwell distribution
+;;; Bonds ;;;
+constraints               = h-bonds      ; controls which bonds in the topology will be converted to rigid holonomic constraints.
+constraint-algorithm      = lincs        ; chooses which solver satisfies any non-SETTLE holonomic constraints. "lincs" is faster but does not work with angle constraints
+lincs-order               = 4            ; accuracy of lincs: the number of matrices in the expansion for the matrix inversion
+lincs-iter                = 1            ; number of iterative corrections to matrix inversion to compensate for lengthening due to rotation
+lincs-warnangle           = 30           ; print warning to log file and stderr if bond rotations be more than this angle
+;;; Vacuum parameters ;;;
+nwall                     = 0            ; number of walls
+wall-type                 = 9-3          ; type of wall potential
+wall-atomtype             = WR WL        ; the atom type name in the force field for each wall.
+wall-density              = 50 50        ; [nm-3] / [nm-2] the number density of the atoms for each wall 

parameters/mdp/molecular-dynamics/nvt_prod_1000ns.mdp

@@ -0,0 +1,60 @@
+; Created by Alec Glisman (GitHub: @alec-glisman) on December 31st, 2021
+;;; GROMACS Documentation: https://manual.gromacs.org/documentation/current/user-guide/mdp-options.html ;;;
+;;; Defines ;;;
+; define                   = -DPOSRES_CRYSTAL ; defines a macro that can be used in the topology file
+;;; Run control ;;;
+integrator                = md           ; "steep" for ener minimization and "md" for a leap frog algorithm for integrating Newton's eq of motion
+dt                        = 0.002        ; 2 femtoseconds
+nsteps                    = 500000000    ; 1000 nanoseconds
+comm-mode                 = Linear       ; "Linear" (remove center for mass translation) "Angular" ( >> >> and rotation around the center of mass) "None" (no restric) 
+nstcomm                   = 100          ; [steps] frequency for center of mass motion removal
+comm-grps                 = Aqueous Crystal   ; groups for which the center of mass motion is removed
+;;; Output control ;;;
+nstxout                   = 0           ; [steps] number of steps that elapse between writing coordinates to the output trajectory file (trr)
+nstvout                   = 0           ; [steps] number of steps that elapse between writing velocities to the output trajectory file (trr)
+nstfout                   = 0           ; [steps] number of steps that elapse between writing forces to the output trajectory file (trr)
+nstcalcenergy             = 100         ; [steps] number of steps that elapse between calculating the energies, 0 is never
+nstenergy                 = 1000        ; [steps] number of steps that elapse between writing energies to energy file
+nstlog                    = 1000        ; [steps] number of steps that elapse between writing energies to the log file
+nstxout-compressed        = 1000        ; [steps] number of steps that elapse between writing to .xtc file
+;;; Neighbor searching ;;;
+cutoff-scheme             = Verlet       ; algorithm to generate neigbor list
+nstlist                   = 10           ; [steps] frequency to update the neighbor list
+pbc                       = xyz          ; periodic boundary conditions
+rlist                     = 1.2000       ; [nm] cut-off distance for the short-range neighbor list
+;;; Electrostatics ;;;
+coulombtype               = PME          ; algorithm to generate electrostatics
+rcoulomb                  = 1.2000       ; [nm] The distance for the Coulomb cut-off
+;;; Van der Waals ;;;
+vdw-type                  = cut-off      ; algorithm to generate Van der Waals
+rvdw                      = 1.2000       ; [nm] distance for the LJ or Buckingham cut-off
+DispCorr                  = AllEnerPres  ; corrections to apply for long-ranged energy and/or pressure
+;;; Ewald ;;;
+fourierspacing            = 0.12         ; [nm] for ordinary Ewald, the ratio of the box dimensions and the spacing determines a lower bound for the number of wave vectors to use in each (signed) direction
+pme-order                 = 4            ; interpolation order for PME. 4 equals cubic interpolation. 
+ewald-rtol                = 1.0e-5       ; the relative strength of the Ewald-shifted direct potential at rcoulomb
+ewald-geometry            = 3d           ; "3d" is the ewald sum which is performed in all the three dimensions or "3dc" the reciprocal sum is still performed in 3D but a forec and potential correction applied in the z dimension to produce a pseudo 2D summation
+;;; Temperature coupling ;;;
+Tcoupl                    = nose-hoover            ; Temp coupling using velocity rescaling. Temperature coupling using a Nose-Hoover extended ensemble. The reference temperature and coupling groups are selected as above, but in this case tau-t controls the period of the temperature fluctuations at equilibrium, which is slightly different from a relaxation time. For NVT simulations the conserved energy quantity is written to the energy and log files.
+tc-grps                   = System                 ; groups to couple separately to Temp bath
+tau-t                     = 0.50                   ; [ps] time constant for coupling (one for each group in tc-grps)
+ref-t                     = 300                    ; [K] reference temperature for coupling (one for each group in tc-grps)
+;;; Pressure Coupling ;;;
+Pcoupl                    = no           ; no pressure coupling in NVT
+;;; Velocity generation ;;;
+gen-vel                   = no           ; "no" do not generate velocities "yes" generate velocities in grompp at temp gen-temp which is only meanigful with integrator "md"
+gen-temp                  = 300          ; [K] temperature for Maxwell distribution
+;;; Bonds ;;;
+constraints               = h-bonds      ; controls which bonds in the topology will be converted to rigid holonomic constraints.
+constraint-algorithm      = lincs        ; chooses which solver satisfies any non-SETTLE holonomic constraints. "lincs" is faster but does not work with angle constraints
+lincs-order               = 4            ; accuracy of lincs: the number of matrices in the expansion for the matrix inversion
+lincs-iter                = 1            ; number of iterative corrections to matrix inversion to compensate for lengthening due to rotation
+lincs-warnangle           = 30           ; print warning to log file and stderr if bond rotations be more than this angle
+;;; Frozen groups ;;;
+freezegrps                = Frozen       ; groups that are to be frozen
+freezedim                 = Y Y Y        ; directions in which the groups are frozen
+;;; Vacuum parameters ;;;
+nwall                     = 0            ; number of walls
+wall-type                 = 9-3          ; type of wall potential
+wall-atomtype             = WR WL        ; the atom type name in the force field for each wall.
+wall-density              = 50 50        ; [nm-3] / [nm-2] the number density of the atoms for each wall

scripts/method/equilibration.sh

@@ -216,9 +216,6 @@ else
             -o "${sim_name}.tpr"
         rm "${previous_sim_name}.gro" || exit 1
 
-        # plumed performance
-        export PLUMED_NUM_THREADS="${CPU_THREADS}"
-
         # run NPT equilibration
         # shellcheck disable=SC2086
         "${GMX_BIN}" -nocopyright mdrun -v \
@@ -388,6 +385,9 @@ else
             sed -i 's/#RESTART/RESTART/g' "plumed.dat" || exit 1
         fi
 
+        # plumed performance
+        export PLUMED_NUM_THREADS="${CPU_THREADS}"
+
         # run production equilibration
         # shellcheck disable=SC2086
         "${GMX_BIN}" -nocopyright mdrun -v \

submission/input/0-goddard-tests/1.0_1PAA_0Ca_104calcite_2nm_crystal_5nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/2.0_8PAA_0Ca_104calcite_5nm_crystal_6nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/3.0_32PAA_32Ca_104calcite_12nm_crystal_12nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/4.0_3PAA_0Ca_104calcite_2nm_crystal_3nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/5.0_3PAA_0Ca_104calcite_1nm_crystal_3nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/5.0_3PAA_0Ca_104calcite_1nm_crystal_4nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/6.0_4PAA_0Ca_104calcite_1nm_crystal_3nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/6.0_4PAA_0Ca_104calcite_1nm_crystal_4nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/0-goddard-tests/7.0_4PAA_0Ca_104calcite_2nm_crystal_4nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='100' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/1-model-system/4.0_16PAA_0Ca_104calcite_9nm_crystal_10nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/1-model-system/4.1_16PAA_0Ca_104calcite_9nm_crystal_10nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/1-model-system/4.2_16PAA_0Ca_104calcite_9nm_crystal_10nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/1-model-system/4.3_16PAA_0Ca_104calcite_9nm_crystal_10nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/1-model-system/4.4_16PAA_0Ca_104calcite_9nm_crystal_10nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acr" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}

submission/input/1-model-system/5.0_16PAAn_0Ca_104calcite_9nm_crystal_10nm_height.sh

@@ -20,7 +20,7 @@ export TEMPERATURE_K='300'       # temperature in Kelvin
 export PRESSURE_BAR='1'          # pressure in bar
 
 # integration
-export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500}
+export INTEGRATION_NS='500' # [ns] final simulation time for production run {100, 500, 1000}
 
 # polyelectrolyte chemistry
 export MONOMER="Acn" # Dominant monomer: {Acr, Acn, Asp, Glu, Ace, Alc}