Hi gromacs users, I'm doing a relative free energy calculation of a ligand-protein complex. I'm transforming a hydrogen atom into methyl. I have generated input topology with FESetup. When I try to run minimization I get that error:
""" Energy minimization has stopped, but the forces have not converged to the Requested precision Fmax It stopped because the algorithm tried to make a new step whose size was too Small, or there was no change in the energy since last step. Either way, we Regard the minimization as converged to within the available machine Precision, given your starting configuration and EM parameters. Double precision normally gives you more accuracy, but this is often not Needed for preparing to run molecular dynamics. You may need to increase your constraint accuracy, or turn Off constraints altogether (set constraints = none in mdp file) Steepest Descents converged to machine precision in 50 steps, But did not reach the requested Fmax <100. Potential Energy = 1.1698845e + 07 Maximum force = 7.4549694e + 05 is atom 5745 Norm of force = 3.5688174e + 03 """ As you can see the Potential Energy is too high and if a look at the system I got is completely messed up. Exactly the same complex with the same parameters but a different perturbation (H -> Cl) gave me no problems. So the minimization crash should be related to the presence of dummy atoms. I read some discussion about the same topic but I was not able to find a solution. Everything looks good to me, I do not see anything wrong in the file .gro or .top. I have tried also to run the same minimization on T4 lysozyme tutorial input (https://ccpforge.cse.rl.ac.uk/gf/project/ccpbiosim/wiki/? pagename=T4+lysozyme), and when there are dummy atom involved I obtain the same problem. Probably I'm missing something in the .mdp or there is some error I can not see. I tried to modify constraint but nothing changed. Any kind of suggestion is really appreciated, thank you in advance. Cheers, Davide Bonanni This is the first part of ligand's topology: """ [ moleculetype ] LIG 3 [ atoms ] ; nr type resno resnm atom cgnr charge mass typeB chargeB massB 1 c 1 LIG C1 1 0.667751 12.0100 c 0.667700 12.0100 2 cc 1 LIG C2 2 -0.469349 12.0100 cc -0.471400 12.0100 3 ca 1 LIG C3 3 0.034751 12.0100 ca 0.039700 12.0100 4 ca 1 LIG C4 4 -0.083249 12.0100 ca -0.032600 12.0100 5 c 1 LIG C5 5 0.742351 12.0100 c 0.742300 12.0100 6 ca 1 LIG C6 6 0.121651 12.0100 ca 0.122600 12.0100 7 ca 1 LIG C7 7 0.029951 12.0100 ca 0.029400 12.0100 8 ca 1 LIG C8 8 0.029951 12.0100 ca 0.029400 12.0100 9 ca 1 LIG C9 9 0.134451 12.0100 ca 0.134400 12.0100 10 ca 1 LIG C10 10 0.134451 12.0100 ca 0.134400 12.0100 11 cp 1 LIG C11 11 -0.137949 12.0100 cp -0.139000 12.0100 12 cp 1 LIG C12 12 -0.012949 12.0100 cp -0.013000 12.0100 13 ca 1 LIG C13 13 -0.108449 12.0100 ca -0.109000 12.0100 14 ca 1 LIG C14 14 -0.108449 12.0100 ca -0.109000 12.0100 15 ca 1 LIG C15 15 -0.136949 12.0100 ca -0.137000 12.0100 16 ca 1 LIG C16 16 -0.136949 12.0100 ca -0.137000 12.0100 17 ca 1 LIG C17 17 -0.134949 12.0100 ca -0.135000 12.0100 18 ca 1 LIG C18 18 -0.155949 12.0100 ca -0.158000 12.0100 19 ca 1 LIG C19 19 -0.122949 12.0100 ca -0.122000 12.0100 20 ca 1 LIG C20 20 -0.226949 12.0100 ca -0.223000 12.0100 21 f 1 LIG F21 21 -0.109349 19.0000 f -0.109400 19.0000 22 f 1 LIG F22 22 -0.130849 19.0000 f -0.131400 19.0000 23 f 1 LIG F23 23 -0.130849 19.0000 f -0.131400 19.0000 24 f 1 LIG F24 24 -0.109349 19.0000 f -0.109400 19.0000 25 nc 1 LIG N25 25 -0.632449 14.0100 nc -0.632500 14.0100 26 na 1 LIG N26 26 0.281751 14.0100 na 0.281700 14.0100 27 n 1 LIG N27 27 -0.465049 14.0100 n -0.464100 14.0100 28 o 1 LIG O28 28 -0.704049 16.0000 o -0.705100 16.0000 29 o 1 LIG O29 29 -0.651049 16.0000 o -0.652100 16.0000 30 h4 1 LIG H30 30 0.158051 1.0080 c3 -0.035800 12.0100 31 ha 1 LIG H31 31 0.141051 1.0080 ha 0.141000 1.0080 32 ha 1 LIG H32 32 0.141051 1.0080 ha 0.141000 1.0080 33 ha 1 LIG H33 33 0.125051 1.0080 ha 0.125000 1.0080 34 ha 1 LIG H34 34 0.125051 1.0080 ha 0.125000 1.0080 35 ha 1 LIG H35 35 0.121051 1.0080 ha 0.121000 1.0080 36 ha 1 LIG H36 36 0.167051 1.0080 ha 0.168000 1.0080 37 ha 1 LIG H37 37 0.117051 1.0080 ha 0.117000 1.0080 38 ha 1 LIG H38 38 0.124051 1.0080 ha 0.124000 1.0080 39 hn 1 LIG H39 39 0.371551 1.0080 hn 0.371500 1.0080 40 du 1 LIG DU40 40 0.000000 1.0000 hc 0.047367 1.0080 41 du 1 LIG DU41 41 0.000000 1.0000 hc 0.047367 1.0080 42 du 1 LIG DU42 42 0.000000 1.0000 hc 0.047367 1.0080 ... ... ... """ This is the em.mdp: """ ; Run control integrator = steep nsteps = 5000 ; EM criteria and other stuff emtol = 100 emstep = 0.01 niter = 20 nbfgscorr = 10 ; Output control nstlog = 1 nstenergy = 1 ; Neighborsearching and short-range nonbonded interactions cutoff-scheme = verlet nstlist = 1 ns_type = grid pbc = xyz rlist = 1.2 ; Electrostatics coulombtype = PME rcoulomb = 1.2 ; van der Waals vdwtype = cutoff vdw-modifier = potential-switch rvdw-switch = 1.0 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; EWALD/PME/PPPM parameters pme_order = 4 ewald_rtol = 1e-5 epsilon_surface = 0 ; Temperature and pressure coupling are off during EM tcoupl = no pcoupl = no ; Free energy control stuff free_energy = yes init_lambda_state = 0 delta_lambda = 0 ; Vectors of lambda specified here ; Each combination is an index that is retrieved from init_lambda_state for each simulation ; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 fep-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 mass-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 ; Options for the decoupling nstdhdl = 100 calc-lambda-neighbors = -1 sc-alpha = 0.5 sc-coul = yes sc-power = 1 sc-r-power = 6 sc-sigma = 0.3 dhdl-derivatives = yes dhdl-print-energy = no separate-dhdl-file = yes dh_hist_size = 0 dh_hist_spacing = 0.1 ; No velocities during EM gen_vel = no ; options for bonds constraints = all-bonds ; Type of constraint algorithm constraint-algorithm = lincs ; Do not constrain the starting configuration continuation = no ; Highest order in the expansion of the constraint coupling matrix lincs-order = 12 lincs-warnangle = 30 ; maximum angle that a bond can rotate energygrps = Protein Non-Protein """ -- 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.
