Hello! I am trying to calculate solvation energy, using Free Energy Perturbation for lutein in water and bilayer. Simulation runs smoothly, but Coulomb part of calculated energy is extremly high for hydrophobic molecule (C40H54(OH)2): above 40kcal/mol // 180 kJ/mol both in water and bilayer. I am working with topology, with OPLS-AA charges.
To check possibility if the charges are somehow wrong I've parametrized molecule with automatic server for charmm - cgenff. Parametrization of the molecule in charmm clearly isn't perfect, but since what I want is only rough estimation of possibility of wrong OPLS-AA charges it should be sufficient. However, I am still getting even higher value of 300kJ/mol. I've also checked methanol under same conditions in mdp file and obtained results were reasonable. Now, I suspect I've messed something in mdp FEP settings,could you please look at it and give me a hint, what may be wrong with these? Best, Krzysiek ; Run control integrator = sd ; Langevin dynamics tinit = 0 dt = 0.002 nsteps = 1750000 ; 3.5 ns nstcomm = 100 ; Output control nstxout = 500 nstvout = 500 nstfout = 0 nstlog = 500 nstenergy = 500 nstxout-compressed = 0 ; Neighborsearching and short-range nonbonded interactions cutoff-scheme = verlet nstlist = 20 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.1 ; EWALD/PME/PPPM parameters pme_order = 5 ewald_rtol = 1e-06 epsilon_surface = 0 ; Temperature coupling ; tcoupl is implicitly handled by the sd integrator tc_grps = System tau_t = 0.6 ref_t = 310 ; Pressure coupling is on for NPT Pcoupl = Parrinello-Rahman pcoupltype = isotropic 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 12 13 14 15 16 17 18 19 20 vdw_lambdas = 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.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 ; We are not transforming any bonded or restrained interactions bonded_lambdas = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 restraint_lambdas = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ; Masses are not changing (particle identities are the same at lambda = 0 and lambda = 1) mass_lambdas = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ; Not doing simulated temperting here temperature_lambdas = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ; Options for the decoupling sc-alpha = 0.5 sc-coul = no ; linear interpolation of Coulomb (none in this case) sc-power = 1 sc-sigma = 0.3 couple-moltype = lutein ; name of moleculetype to decouple couple-lambda0 = vdw-q ; all couple-lambda1 = vdw ; vdw couple-intramol = yes nstdhdl = 100 ; No velocities during EM gen_vel = no ; options for bonds constraints = h-bonds ; we only have C-H bonds here ; 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 = 4 -- Jagiellonian University Department of Computational Biophysics and Bioinformatics tel.1: (12) 664 61 49 tel.2: (48) 664 086 049 -- 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.