GROMACS version:
GROMACS modification: Yes/No
Here post your question
Hello again,
I have a rather large system of 60 proteins and I know where the overlapping is coming from however it is a part of the system. How do I change to soft-core potential with the free energy code? I am new to command line and coding, so the manual doesn’t quite make sense to me about what I should change in the either the mdp file to run soft core potentials or if I should add a command to execute this to see if my system can be minimized. Can you explain to me? Thanks
Heather
Command line:
gmx mdrun -v -deffnm em
Reading file em.tpr, VERSION 2023 (single precision)
Using 64 MPI threads
Using 1 OpenMP thread per tMPI thread
Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = 50000
Step= 0, Dmax= 1.0e-02 nm, Epot= 7.07618e+19 Fmax= inf, atom= 139826
Step= 14, Dmax= 1.2e-06 nm, Epot= 7.07618e+19 Fmax= inf, atom= 139826
Energy minimization has stopped because the force on at least one atom is not
finite. This usually means atoms are overlapping. Modify the input
coordinates to remove atom overlap or use soft-core potentials with the free
energy code to avoid infinite forces.
You could also be lucky that switching to double precision is sufficient to
obtain finite forces.
writing lowest energy coordinates.
Steepest Descents converged to machine precision in 15 steps,
but did not reach the requested Fmax < 1000.
Potential Energy = 7.0761820e+19
Maximum force = inf on atom 139826
Norm of force = inf
Hi,
I think you should first try to tune the parameter in the mdp file for energy minimization.
Or you can try to rebuild the starting structure, or try to understand if the issue is in the structure of one protein or the interaction between the proteins
\Alessandra
Hello,
I changed the n-step to 100, it still ran into the error. I think the problem is the structure as it is a capsid. They are all the same viral protein with an icosahedral structure, so I think the interaction is between the proteins. If I were to try the free energy, would I add it to the mdp file under parameters with the free energy = yes ?
Heather
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 100 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
free-energy = yes ; soft core potentials used
Like this?
The free-energy keyword activates the free energy code, but you need to specify what entity you want to scale, whether you want vdW and/or electrostatics scaled, and by how much. Use of the free energy code to relax bad contacts is a rather complicated way to resolve clashes, when in most cases a simple in vacuo minimization works quite well.
If you want to pursue the free energy option, the basics of the relevant keywords are discussed here: Free Energy Calculations You would want to scale the entire system and only use one λ state rather than a complete transformation.