bug fixes: two-slab system functional

alec-glisman · Apr 23, 2024 · 165b786 · 165b786
1 parent 9f14a51
commit 165b786
Show file tree

Hide file tree

Showing 8 changed files with 64 additions and 68 deletions.
diff --git a/initial-structure/calcium-carbonate-crystal/generation/python/requirements.yml b/initial-structure/calcium-carbonate-crystal/generation/python/requirements.yml
@@ -29,6 +29,7 @@ dependencies:
   - pymatgen
   - openbabel
   - rdkit
+  - gromacswrapper
   # plotting
   - matplotlib
   - seaborn

diff --git a/parameters/mdp/energy-minimization/em.mdp b/parameters/mdp/energy-minimization/em.mdp
@@ -1,15 +1,15 @@
 ; Created by Alec Glisman (GitHub: @alec-glisman) on March 7th, 2022
 ;;; 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
+; define                     = -DPOSRES_CRYSTAL  ; defines a macro that can be used in the topology file
 ;;; Run control ;;;
 integrator                = steep             ; "steep" for energy minimization and "md" for a leap frog algorithm for integrating Newton's eq of motion
 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
 ;;; Energy minimization ;;;
 emtol                     = 10.0         ; [kJ mol^{-1} nm^{-1}] the minimization is converged when the maximum force is smaller than this value
-emstep                    = 0.20         ; [nm] initial step-size
+emstep                    = 0.10         ; [nm] initial step-size
 nsteps                    = 100000       ; [steps] Maximum number of (minimization) steps to perform
 ;;; Output control ;;;
 nstxout                   = 0            ; [steps] number of steps that elapse between writing coordinates to the output trajectory file (trr)
@@ -18,10 +18,10 @@ nstfout                   = 0            ; [steps] number of steps that elapse b
 nstlog                    = 100          ; [steps] number of steps that elapse between writing energies to the log file
 nstcalcenergy             = 100          ; [steps] number of steps that elapse between calculating the energies, 0 is never
 nstenergy                 = 100          ; [steps] number of steps that elapse between writing energies to energy file
-nstxout-compressed        = 100          ; [steps] number of steps that elapse between writing to .xtc file
+nstxout-compressed        = 10           ; [steps] number of steps that elapse between writing to .xtc file
 ;;; Neighbor searching ;;;
 cutoff-scheme             = Verlet       ; algorithm to generate neigbor list
-nstlist                   = 20           ; [steps] frequency to update the neighbor 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 ;;;
@@ -48,8 +48,8 @@ lincs-order               = 4            ; accuracy of lincs: the number of matr
 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                = Crystal       ; groups that are to be frozen
-freezedim                 = Y Y Y         ; directions in which the groups are frozen
+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

diff --git a/parameters/mdp/molecular-dynamics/npt_eqbm.mdp b/parameters/mdp/molecular-dynamics/npt_eqbm.mdp
@@ -40,7 +40,7 @@ tc-grps                   = System                 ; groups to couple separately
 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.
+Pcoupl                    = C-rescale             ; 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

diff --git a/scripts/method/equilibration.sh b/scripts/method/equilibration.sh
@@ -97,13 +97,6 @@ else
             sed -i 's/^rcoulomb.*/rcoulomb = 0.7/g' "${sim_name}.mdp" || exit 1
             sed -i 's/^rvdw.*/rvdw = 0.7/g' "${sim_name}.mdp" || exit 1
         fi
-        # add vacuum parameters to mdp file
-        if [[ "${VACUUM}" == 'True' ]]; then
-            sed -i 's/^ewald-geometry .*/ewald-geometry            = 3dc/g' "${sim_name}.mdp" || exit 1
-            sed -i 's/^pbc .*/pbc                       = xy/g' "${sim_name}.mdp" || exit 1
-            sed -i 's/^nwall .*/nwall                     = 2/g' "${sim_name}.mdp" || exit 1
-            echo "wall-r-linpot = 0.1" >>"${sim_name}.mdp" || exit 1
-        fi
 
         # make tpr file for NVT equilibration
         "${GMX_BIN}" -nocopyright grompp \
@@ -199,12 +192,6 @@ else
             sed -i 's/^rcoulomb.*/rcoulomb = 0.7/g' "${sim_name}.mdp" || exit 1
             sed -i 's/^rvdw.*/rvdw = 0.7/g' "${sim_name}.mdp" || exit 1
         fi
-        # add vacuum parameters to mdp file
-        if [[ "${VACUUM}" == 'True' ]]; then
-            sed -i 's/^ewald-geometry .*/ewald-geometry            = 3dc/g' "${sim_name}.mdp" || exit 1
-            sed -i 's/^pbc .*/pbc                       = xy/g' "${sim_name}.mdp" || exit 1
-            sed -i 's/^nwall .*/nwall                     = 2/g' "${sim_name}.mdp" || exit 1
-        fi
 
         # make tpr file
         "${GMX_BIN}" -nocopyright grompp \
@@ -247,6 +234,15 @@ EOF
             --gro_filename "${sim_name}.gro" \
             --percentile '0.4'
 
+        # unwrap final gro file
+        "${GMX_BIN}" -nocopyright trjconv \
+            -f "${sim_name}.gro" \
+            -s "${sim_name}.tpr" \
+            -o "${sim_name}.gro" \
+            -pbc 'mol' -ur 'tric' <<EOF
+System
+EOF
+
         # convert final gro file to pdb file
         "${GMX_BIN}" -nocopyright trjconv \
             -f "${sim_name}.gro" \

diff --git a/scripts/method/initialization.sh b/scripts/method/initialization.sh
@@ -106,12 +106,6 @@ echo "INFO: Copying input files to working directory"
         sed -i 's/^rcoulomb.*/rcoulomb = 0.7/g' "mdin.mdp"
         sed -i 's/^rvdw.*/rvdw = 0.7/g' "mdin.mdp"
     fi
-    # add vacuum parameters to mdp file
-    if [[ "${VACUUM}" == "True" ]]; then
-        sed -i 's/^ewald-geometry .*/ewald-geometry            = 3dc/g' "mdin.mdp"
-        sed -i 's/^pbc .*/pbc                       = xy/g' "mdin.mdp"
-        sed -i 's/^nwall .*/nwall                     = 2/g' "mdin.mdp"
-    fi
 } >>"${log_file}" 2>&1
 
 # ##############################################################################
@@ -182,8 +176,10 @@ echo "INFO: Importing structure to Gromacs"
     carbonate_carbon_z="$(grep 'CRB' "crystal.gro" | grep -o '.\{6\}$' | awk '{$1=$1};1')"
     minimum_z_coord="$(echo "${carbonate_carbon_z}" | sort -n)"
     z_min="$(echo "${minimum_z_coord}" | awk 'NR==1{print $1}')"
+    z_max="$(echo "${minimum_z_coord}" | awk 'END{print $1}')"
+
     # subtract offset [nm] to z_min to ensure that all atoms are within the box and we can see water structure
-    offset='0.0'
+    offset='0.1'
     z_min="$(bc <<<"scale=5; ${z_min} - ${offset}")"
     echo "DEBUG: Minimum z-coordinate of crystal [nm]: ${z_min}"
 
@@ -194,12 +190,14 @@ echo "INFO: Importing structure to Gromacs"
             -f "${sim_name}.gro" \
             -o "${sim_name}.gro" \
             -translate '0' '0' "${z_min}"
+        z_max="$(bc <<<"scale=5; ${z_max} + ${z_min}")"
     else
         # shift z-coordinates of all atoms by z_min
         "${GMX_BIN}" -nocopyright editconf \
             -f "${sim_name}.gro" \
             -o "${sim_name}.gro" \
             -translate '0' '0' "-${z_min}"
+        z_max="$(bc <<<"scale=5; ${z_max} - ${z_min}")"
     fi
 
     # wrap all atoms into box
@@ -230,15 +228,15 @@ echo "INFO: Adding solvent"
 
     # subtract z_min from pdb bulk z-coordinates and add buffer for outer layers
     buffer='0.3'
-    # gro_zmin="$(bc <<<"scale=5; ${PDB_BULK_ZMIN} - ${z_min} - ${buffer}")"
-    gro_zmax="$(bc <<<"scale=5; ${PDB_BULK_ZMAX} - ${z_min} + ${buffer}")"
+    gro_zmin="$(bc <<<"scale=5; (${PDB_BULK_ZMIN} - ${z_min} - ${buffer})")"
+    gro_zmax="$(bc <<<"scale=5; (${PDB_BULK_ZMAX} - ${z_min} + ${buffer})")"
 
     # find "bad" water molecules that are inside the crystal
     "${GMX_BIN}" -nocopyright select \
         -f "${sim_name}.gro" \
         -s "${sim_name}.gro" \
         -on "bad_waters.ndx" <<EOF
-"Bad_SOL" same residue as (name OW HW1 HW2 and (z <= ${gro_zmax}))
+"Bad_SOL" same residue as (name OW HW1 HW2 and (z > ${gro_zmin} and z < ${gro_zmax}))
 EOF
     # remove "f0_t0.000" from index file groups
     sed -i 's/_f0_t0.000//g' "bad_waters.ndx"
@@ -264,12 +262,10 @@ Good_Atoms
 EOF
 } >>"${log_file}" 2>&1
 
-# echo "DEBUG: Minimum z-coordinate of crystal after solvation [nm]: ${gro_zmin}"
-echo "DEBUG: Maximum z-coordinate of crystal after solvation [nm]: ${gro_zmax}"
-
 # find number of "bad" water molecules from log file
 n_bad_atoms="$(grep "Bad_SOL" "${log_file}" | head -n 1 | awk '{print $4}')"
 n_bad_waters="$((n_bad_atoms / 3))"
+echo "DEBUG: Water exclusion zone [nm]: ${gro_zmin} to ${gro_zmax}"
 echo "DEBUG: Number of water molecules inside crystal (removed): ${n_bad_waters}"
 
 {
@@ -517,7 +513,7 @@ echo "INFO: Running energy minimization"
         -f "${sim_name}.trr" \
         -s "${sim_name}.tpr" \
         -o "${sim_name}.gro" \
-        -dump '100000' <<EOF
+        -dump '100000000' <<EOF
 System
 EOF
 
@@ -527,7 +523,7 @@ EOF
         -s "${sim_name}.tpr" \
         -o "${sim_name}_final.pdb" \
         -pbc 'mol' -ur 'tric' \
-        -dump '100000' -conect <<EOF
+        -dump '100000000' -conect <<EOF
 System
 System
 EOF

diff --git a/scripts/utilities/second_slab.sh b/scripts/utilities/second_slab.sh
@@ -9,7 +9,7 @@
 #             submission/input/*.sh script. Script should only be called from
 #             the equilibration.sh script.
 
-z_buffer="0.0" # [nm]
+z_buffer="0.1" # [nm]
 sim_name="npt_eqbm"
 archive_dir="3-second-slab"
 log_file="second_slab.log"
@@ -41,7 +41,7 @@ else
         z_max="$(echo "${z_coords}" | tail -n 1)"
         z_min="$(echo "${z_coords}" | head -n 1)"
         z_height_crystal="$(echo "scale=5; ${z_max} - ${z_min}" | bc)"
-        z_translate="$(echo "scale=5; ${z_height} + ${z_max}" | bc)"
+        z_translate="$(echo "scale=5; ${z_height} + ${z_max} + ${z_buffer}" | bc)"
         z_height_2slab="$(echo "scale=5; ${z_height} + ${z_height_crystal} + ${z_buffer}" | bc)"
 
         # insert "-" at column 39 in all lines of gro file to reflect crystal on z-axis

diff --git a/...1.0_104calcite_7nm_crystal_10nm_height.sh → .../6.0_104calcite_9nm_crystal_4nm_height.sh b/...1.0_104calcite_7nm_crystal_10nm_height.sh → .../6.0_104calcite_9nm_crystal_4nm_height.sh
@@ -15,7 +15,7 @@ export GPU_IDS='-1'     # GPU device(s) to use (0 = first GPU, 01 = first two GP
 # System components ###########################################################
 
 # tag for system
-export TAG_JOBID="8.1.0-00idx" # tag to append to system name
+export TAG_JOBID="7.6.0-00idx" # tag to append to system name
 
 # statistical mechanics
 export PRODUCTION_ENSEMBLE='NVT' # {NVT, NPT}
@@ -35,10 +35,10 @@ export SURFACE="104"     # Miller index of crystal surface {104, 001}
 
 # system size
 export VACUUM='True'        # {True, False}
-export SURFACE_SIZE='7'     # size of crystal surface in nm {3, 5, 8, 9, 10, 11, 12, 13}
-export BOX_HEIGHT='10'      # height of simulation box in nm
-export PDB_BULK_ZMIN='0.46' # z-coordinate of bottom of bulk part of crystal in nm in PDB file
-export PDB_BULK_ZMAX='0.30' # z-coordinate of top of bulk part of crystal in nm in PDB file
+export SURFACE_SIZE='9'     # size of crystal surface in nm {3, 5, 8, 9, 10, 11, 12, 13}
+export BOX_HEIGHT='4'       # height of simulation box in nm
+export PDB_BULK_ZMIN='4.62' # z-coordinate of bottom of bulk part of crystal in nm in PDB file
+export PDB_BULK_ZMAX='5.09' # z-coordinate of top of bulk part of crystal in nm in PDB file
 
 # number of each component
 export N_SLAB='2'      # number of crystal slabs {1, 2}
@@ -79,8 +79,8 @@ export METAD_GRID_MAX='10'        # maximum grid point for Gaussian deposition
 export METAD_PACE='500'           # [steps] between deposition of Gaussians
 
 # harmonic restraints
-export ATOM_REFERENCE='1708'  # atom number of reference atom for harmonic restraints (1 = first atom)
-export ATOM_OFFSET='-0.305'   # z-coordinate offset of reference atom from crystal surface in nm
 export PE_WALL_MIN='0.3'      # z-coordinate of lower wall in nm
 export PE_WALL_MAX='4.0'      # z-coordinate of upper wall in nm
 export PE_WALL_MAX_EQBM='2.0' # z-coordinate of upper wall in nm during equilibration
+export ATOM_REFERENCE='7310'  # atom number of reference atom for harmonic restraints (1 = first atom)
+export ATOM_OFFSET='-0.305'   # z-coordinate offset of reference atom from crystal surface in nm
diff --git a/..._110aragonite_12nm_crystal_10nm_height.sh → ...2.1_010aragonite_16paa_0Ca_12nmcrystal.sh b/..._110aragonite_12nm_crystal_10nm_height.sh → ...2.1_010aragonite_16paa_0Ca_12nmcrystal.sh
@@ -15,7 +15,7 @@ export GPU_IDS='-1'     # GPU device(s) to use (0 = first GPU, 01 = first two GP
 # System components ###########################################################
 
 # tag for system
-export TAG_JOBID="8.2.0-00idx" # tag to append to system name
+export TAG_JOBID="8.2.1" # tag to append to system name
 
 # statistical mechanics
 export PRODUCTION_ENSEMBLE='NVT' # {NVT, NPT}
@@ -31,26 +31,36 @@ export BLOCK=""      # Block copolymer: {iiia, iiiiiiia, iiic, iiiiiic} for {Acr
 
 # calcium carbonate crystal surface
 export CRYSTAL="aragonite" # {calcite, aragonite, vaterite}
-export SURFACE="110"       # Miller index of crystal surface {104, 001}
+export SURFACE="010"       # Miller index of crystal surface {104, 001}
+
+# TODO: update the following parameters
+# calcium carbonate crystalline references
+export PDB_BULK_ZMIN='4.62'  # z-coordinate of bottom of bulk part of crystal in nm in PDB file
+export PDB_BULK_ZMAX='5.09'  # z-coordinate of top of bulk part of crystal in nm in PDB file
+export ATOM_REFERENCE='7310' # atom number of reference atom for harmonic restraints (1 = first atom)
+export ATOM_OFFSET='-0.305'  # z-coordinate offset of reference atom from crystal surface in nm
 
 # system size
-export VACUUM='True'        # {True, False}
-export SURFACE_SIZE='12'    # size of crystal surface in nm {3, 5, 8, 9, 10, 11, 12, 13}
-export BOX_HEIGHT='10'      # height of simulation box in nm
-export PDB_BULK_ZMIN='0.46' # z-coordinate of bottom of bulk part of crystal in nm in PDB file
-export PDB_BULK_ZMAX='0.30' # z-coordinate of top of bulk part of crystal in nm in PDB file
+export VACUUM='True'     # {True, False}
+export SURFACE_SIZE='12' # size of crystal surface in nm {3, 5, 8, 9, 10, 11, 12, 13}
+export BOX_HEIGHT='14'   # height of simulation box in nm
 
 # number of each component
-export N_SLAB='2'      # number of crystal slabs {1, 2}
-export N_MONOMER='0'   # number of monomers in chain {1, 2, 5, 8, 16, 32}
-export N_CHAIN='0'     # number of chains
+export N_SLAB='1'      # number of crystal slabs {1, 2}
+export N_MONOMER='16'  # number of monomers in chain {1, 2, 5, 8, 16, 32}
+export N_CHAIN='1'     # number of chains
 export N_CARBONATE='0' # number of aqueous carbonate ions
-export N_SODIUM='0'    # number of aqueous sodium ions
+export N_SODIUM='16'   # number of aqueous sodium ions
 export N_CALCIUM='0'   # number of aqueous calcium ions
 export N_CHLORINE='0'  # number of aqueous chlorine ions
 
 # Enhanced sampling ###########################################################
 
+# harmonic restraints
+export PE_WALL_MAX_EQBM='2.0' # z-coordinate of upper wall in nm during equilibration
+export PE_WALL_MIN='0.3'      # z-coordinate of lower wall in nm
+export PE_WALL_MAX='4.0'      # z-coordinate of upper wall in nm
+
 # hamiltonian replica exchange
 export HREMD_N_REPLICA='32' # number of replicas in HREMD simulations
 export HREMD_T_MIN='300'    # [K] min temperature of HREMD replicas, set to TEMPERATURE_K
@@ -61,12 +71,12 @@ export HREMD_N_STEPS='1000' # number of steps between replica exchange attempts
 export ONEOPES_N_REPLICA='8'        # number of replicas in OneOPES simulations
 export ONEOPES_N_STEPS='1000'       # number of steps between replica exchange attempts
 export ONEOPES_LARGE_BARRIER='30'   # [kJ/mol] large barrier height for OneOPES replica exchange
-export ONEOPES_SMALL_BARRIER='5'    # [kJ/mol] small barrier height for OneOPES replica exchange
-export ONEOPES_REPLICA_2_TEMP='300' # [K] max OPES MultiTherm temperature of replica 2
-export ONEOPES_REPLICA_3_TEMP='310' # [K] max OPES MultiTherm temperature of replica 3
-export ONEOPES_REPLICA_4_TEMP='325' # [K] max OPES MultiTherm temperature of replica 4
-export ONEOPES_REPLICA_5_TEMP='340' # [K] max OPES MultiTherm temperature of replica 5
-export ONEOPES_REPLICA_6_TEMP='355' # [K] max OPES MultiTherm temperature of replica 6
+export ONEOPES_SMALL_BARRIER='8'    # [kJ/mol] small barrier height for OneOPES replica exchange
+export ONEOPES_REPLICA_2_TEMP='304' # [K] max OPES MultiTherm temperature of replica 2
+export ONEOPES_REPLICA_3_TEMP='312' # [K] max OPES MultiTherm temperature of replica 3
+export ONEOPES_REPLICA_4_TEMP='326' # [K] max OPES MultiTherm temperature of replica 4
+export ONEOPES_REPLICA_5_TEMP='338' # [K] max OPES MultiTherm temperature of replica 5
+export ONEOPES_REPLICA_6_TEMP='353' # [K] max OPES MultiTherm temperature of replica 6
 export ONEOPES_REPLICA_7_TEMP='370' # [K] max OPES MultiTherm temperature of replica 7
 
 # well-tempered metadynamics
@@ -77,10 +87,3 @@ export METAD_GRID_SPACING='0.005' # width of bins in the meta-dynamics grid
 export METAD_GRID_MIN='0'         # minimum grid point for Gaussian deposition
 export METAD_GRID_MAX='10'        # maximum grid point for Gaussian deposition
 export METAD_PACE='500'           # [steps] between deposition of Gaussians
-
-# harmonic restraints
-export ATOM_REFERENCE='1708'  # atom number of reference atom for harmonic restraints (1 = first atom)
-export ATOM_OFFSET='-0.305'   # z-coordinate offset of reference atom from crystal surface in nm
-export PE_WALL_MIN='0.3'      # z-coordinate of lower wall in nm
-export PE_WALL_MAX='4.0'      # z-coordinate of upper wall in nm
-export PE_WALL_MAX_EQBM='2.0' # z-coordinate of upper wall in nm during equilibration