Hi, everyone, I am running into an issue where grompp is using a tremendous amount of memory and crashing, even though my system is not especially large (63976 atoms).
I am using GROMACS 2016.3. My system consists of liquid water (7,930 molecules) next to a block of ice (8,094 molecules). The ice oxygens are restrained to their lattice position with a harmonic potential with strength k = 4,000 kJ/mol/nm^2. I am using the TIP4P/Ice model, which is a rigid 4-site model with a negative partial charge located on a virtual site rather than the oxygen. My goal is to systematically reduce the electrostatic interactions between the water molecules and the position-restrained ice, while leaving water-water and ice-ice interactions unaffected. To accomplish this, I am modeling all of the ice molecules using a single moleculetype so that I can take advantages of GROMACS' FEP features to selectively scale interactions. I explicitly specify all constraints and exclusions in the topology file. This moleculetype contains one virtual site, 3 constraints, and 4 exclusions per "residue" (ice molecule). When I run grompp, I receive the following error, which I think means that a huge block of memory (~9 GB) was requested but could not be allocated: ===== Command line: gmx grompp -f npt.mdp -c md.gro -p topol.top -n index.ndx -r initconf_packmol.gro -o input.tpr -maxwarn 2 -pp processed.top ... Generated 21 of the 21 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 21 of the 21 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'ICE' turning H bonds into constraints... Excluding 3 bonded neighbours molecule type 'SOL' turning H bonds into constraints... Coupling 1 copies of molecule type 'ICE' Setting gen_seed to 1021640799 Velocities were taken from a Maxwell distribution at 273 K Cleaning up constraints and constant bonded interactions with virtual sites Removing all charge groups because cutoff-scheme=Verlet ------------------------------------------------------- Program: gmx grompp, version 2016.3 Source file: src/gromacs/utility/smalloc.cpp (line 226) Fatal error: Not enough memory. Failed to realloc -8589934588 bytes for il->iatoms, il->iatoms=25e55010 (called from file /home/semarks/source/gromacs/2016.3/icc/serial/gromacs-2016.3/src/gromacs/ gmxpreprocess/convparm.cpp, line 565) For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors ------------------------------------------------------- ======= In the hope that it helps with diagnosing the problem, here is my mdp file: ======= ; Water and ice ; RUN CONTROL PARAMETERS ; Molecular dynamics integrator = md ; Leapfrog dt = 0.002 ; 2 fs - time step [ps] nsteps = 5000000 ; 10 ns - run time [steps] ; Center of mass motion (COMM) removal comm-mode = none ; Linear nstcomm = 10 ; Removal frequency (>= nstcalcenergy) [steps] comm-grps = System ; Groups for COMM removal (blank => whole system) ; Initial velocities gen_vel = yes ; Generate velocities using Boltzmann distribution gen_temp = 273 ; Temperature for Boltzmann distribution [K] gen_seed = -1 ; Seed for RNG is the job ID ; OUTPUT CONTROL OPTIONS nstcalcenergy = 1 ; Frequency of energy calculation [steps] ; Output frequency for coords (x), velocities (v) and forces (f) ; trr nstxout = 0 ; never - print coordinates [steps] nstvout = 0 ; never - print velocities [steps] nstfout = 0 ; never - print forces [steps] ; Log end edr files nstlog = 2500 ; 5 ps - print energies to log file [steps] nstenergy = 2500 ; 5 ps - print energies to energy file [steps] ; xtc nstxout-compressed = 2500 ; 5 ps - print coordinates to xtc file [steps] compressed-x-precision = 1000 ; Number of zeros = number of places after decimal point compressed-x-grps = System ; Groups to write to xtc file ; BOUNDARY CONDITIONS pbc = xyz periodic_molecules = no ; Rigid graphene has no intramolecule interactions ; NEIGHBOR LIST cutoff-scheme = Verlet nstlist = 10 ; Neighbor list update frequency [steps] ns-type = grid ; More efficient than simple search verlet-buffer-tolerance = 0.005 ; [kJ/mol/ps] (-1.0 --> use rlist) rlist = 1.0 ; Cutoff distance for short-range neighbor list [nm] ; VAN DER WAALS vdwtype = Cut-off vdw-modifier = Potential-shift-Verlet rvdw = 1.0 ; Radius for vdW cutoff [nm] DispCorr = no ; ELECTROSTATICS coulombtype = PME ; Fast, smooth, particle mesh Ewald (PME) coulomb-modifier = Potential-shift-Verlet rcoulomb = 1.0 ; Short-range Coulomb cutoff [nm] epsilon-r = 1 ; Relative dielectric constant ; Ewald/PME/PPPM parameters fourierspacing = 0.12 ; FFT grid spacing [nm] pme_order = 4 ; Cubic interpolation of PME grid ewald-rtol = 1e-05 ; ewald-geometry = 3d ; 3-D charge treatment epsilon-surface = 0 ; No dipole correction ; THERMOSTAT Tcoupl = V-rescale ; stochastic velocity-rescale nstTcouple = 1 tc-grps = ICE SOL ; Groups which are controlled together tau_T = 0.5 0.5 ; Time constant for coupling [ps] ref_T = 273.0 273.0 ; Temperature setpoint [K] ; BAROSTAT Pcoupl = Parrinello-Rahman ; (options: no/Berendsen/Parrinello-Rahman) Pcoupltype = anisotropic nstPcouple = 1 tau-P = 4.0 ; Time constant for coupling [ps] ref-P = 1.0 1.0 1.0 0.0 0.0 0.0 ; Pressure setpoint [bar] compressibility = 0.0 4.8e-5 0.0 0.0 0.0 0.0 ; Approximate compressibility [bar^-1] refcoord-scaling = no ; CONSTRAINTS constraints = H-bonds ; Convert bonds with H to constraints constraint_algorithm = LINCS ; Best choice if you don't have angle constraints lincs-order = 4 ; Highest order in expansion of constraint coupling matrix lincs-iter = 1 ; More iterations --> higher accuracy continuation = no ; ; Use GROMACS' FEP features to scale ice-water electrostatics, while leaving ; ice-ice electrostatics unaffected. ; - This setup: lambda = scaling factor for interactions (ranges from 0 to 1) free-energy = yes ; turns on FEP couple-moltype = ICE ; group for coupling couple-intramol = no ; don't adjust ice-ice interactions couple-lambda0 = vdw ; ice-water electrostatics are off at lambda=0 couple-lambda1 = vdw-q ; all interactions are on at lambda=1 init-lambda = 0.5 ; value of lambda (stays fixed with this setup) delta-lambda = 0.0 ; don't change lambda nstdhdl = 0 ; no extra output needed dhdl-derivatives = no ; no extra output needed ======= Thank you for your time. Sincerely, Sean -- Sean M. Marks PhD Candidate Dept. of Chemical and Biomolecular Engineering University of Pennsylvania -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.