From 165b786d14bfb3ce536b2cd0a4419f8fa019eddd Mon Sep 17 00:00:00 2001 From: Alec <30010253+alec-glisman@users.noreply.github.com> Date: Tue, 23 Apr 2024 15:34:26 -0700 Subject: [PATCH] bug fixes: two-slab system functional --- .../generation/python/requirements.yml | 1 + parameters/mdp/energy-minimization/em.mdp | 12 ++--- .../mdp/molecular-dynamics/npt_eqbm.mdp | 2 +- scripts/method/equilibration.sh | 22 ++++---- scripts/method/initialization.sh | 26 ++++------ scripts/utilities/second_slab.sh | 4 +- .../6.0_104calcite_9nm_crystal_4nm_height.sh} | 14 ++--- ...2.1_010aragonite_16paa_0Ca_12nmcrystal.sh} | 51 ++++++++++--------- 8 files changed, 64 insertions(+), 68 deletions(-) rename submission/input/{8-two-slab/1.0_104calcite_7nm_crystal_10nm_height.sh => 7-vacuum-systems/6.0_104calcite_9nm_crystal_4nm_height.sh} (91%) rename submission/input/8-two-slab/{2.0_110aragonite_12nm_crystal_10nm_height.sh => 2.1_010aragonite_16paa_0Ca_12nmcrystal.sh} (70%) diff --git a/initial-structure/calcium-carbonate-crystal/generation/python/requirements.yml b/initial-structure/calcium-carbonate-crystal/generation/python/requirements.yml index 679e631..bab882a 100644 --- 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 index f4ca819..fb8747d 100644 --- a/parameters/mdp/energy-minimization/em.mdp +++ b/parameters/mdp/energy-minimization/em.mdp @@ -1,7 +1,7 @@ ; 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) @@ -9,7 +9,7 @@ nstcomm = 100 ; [steps] frequency for center of 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 index 8871727..9146cbd 100644 --- 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 index e19342f..0e24aed 100755 --- 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' <>"${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" < ${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' <