Dear Gromacs users, I am using Gromacs 5.0.6 with CHARMM27 forcefield. I want to do FEP analysis and I did it based on this tutorials: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/01_theory.html
. But I just want to peturb some atoms of the protein, however, in the mdp file,I only can choose whole protein moleculetype. couple-moltype = Protein I don’t know how to calculate free energy by only perturbing charge of some atoms of protein. and I have tried to modify topology file that I added the state B, liking the following: but it seems unable to recognize the state B in topology file. 25 C 3 ALA C 25 0.51 12.011 C 0.00 12.011 ; qtot 0.465 26 O 3 ALA O 26 -0.51 15.999 O 0.00 15.999 ; qtot -0.045 ; residue 4 ALA rtp ALA q 0.0 27 NH1 4 ALA N 27 -0.47 14.007 NH1 0.00 14.007 ; qtot -0.515 28 H 4 ALA HN 28 0.31 1.008 H 0.00 1.008 ; qtot -0.205 29 CT1 4 ALA CA 29 0.07 12.011 CT1 0.00 12.011 ; qtot -0.135 30 HB 4 ALA HA 30 0.09 1.008 HB 0.00 1.008 ; qtot -0.045 31 CT3 4 ALA CB 31 -0.27 12.011 ; qtot -0.315 32 HA 4 ALA HB1 32 0.09 1.008 ; qtot -0.225 33 HA 4 ALA HB2 33 0.09 1.008 ; qtot -0.135 34 HA 4 ALA HB3 34 0.09 1.008 ; qtot -0.045 35 C 4 ALA C 35 0.51 12.011 ; qtot 0.465 36 O 4 ALA O 36 -0.51 15.999 ; qtot -0.045 ; residue 5 ALA rtp ALA q 0.0 37 NH1 5 ALA N 37 -0.47 14.007 ; qtot -0.515 38 H 5 ALA HN 38 0.31 1.008 ; qtot -0.205 but it seems that it can't recognize state B in topology file. And it seems to calculate the free energy just based on the mdp file, couple-moltype = Protein Mdp file is the following: ; Pressure coupling is on for NPT Pcoupl = Parrinello-Rahman tau_p = 1.0 compressibility = 4.5e-05 ref_p = 1.0 ; Free energy control stuff free_energy = yes init_lambda_state = 0 delta_lambda = 0 calc_lambda_neighbors = 1 ; only immediate neighboring windows ; 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 vdw_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 coul_lambdas = 0.0 0.01 0.02 0.04 0.06 0.08 0.1 0.15 0.25 0.5 0.75 1.0 ; We are not transforming any bonded or restrained interactions bonded_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 restraint_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 ; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1) 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 ; Not doing simulated temperting here temperature_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 ; Options for the decoupling sc-alpha = 0.5 sc-coul = no ; linear interpolation of Coulomb (none in this case) sc-power = 1.0 sc-sigma = 0.3 couple-moltype = Protein ; name of moleculetype to decouple couple-lambda0 = vdw-q ; only van der Waals interactions couple-lambda1 = none ; turn off everything, in this case only vdW couple-intramol = no nstdhdl = 10 ; Do not generate velocities gen_vel = no ; options for bonds constraints = h-bonds ; we only have C-H bonds here ; Type of constraint algorithm constraint-algorithm = lincs ; Constrain the starting configuration ; since we are continuing from NPT continuation = yes ; Highest order in the expansion of the constraint coupling matrix lincs-order = 12 I’m wandering how can I solve the problem ? It would be a great help if anyone could help me out here. Many thanks! Zhang Xiao
-- 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.
