Re: MMGBSA binding free energy by MolAICal and NAMD

From: molaical (molaical_at_aliyun.com)
Date: Tue Jul 27 2021 - 12:36:12 CDT

Dear friend,

Firstly, I think you should check whether your receptor and ligand are always in the same image along simulated time. The ligand maybe have run out of this periodic mirror image. Please check this problem. If the ligand is out of the current phase, please use the PBC of VMD to wrap the ligand into receptor in the current phase.

Thanks.

------------------------------------------------------------------
Sender:Mi Yang <drmiyang2019_at_gmail.com>
Sent At:2021 Jul. 15 (Thu.) 02:44
Recipient:"<namd-l_at_ks.uiuc.edu>" <namd-l_at_ks.uiuc.edu>
Subject:namd-l: MMGBSA binding free energy by MolAICal and NAMD

Dear Colleagues,
We have calculated MM/GBSA energy using MolAICal and NAMD. Its final value is coming positive even for nicely bounded poses after MD. Any idea what parameters can be adjusted in MD inputs. All inputs are given below:

=======MolAICal out with MM/GBSA energy=======
Total run 1010 frames from complex.log.
Average BOND: 952.4267987128711
Average ANGLE: 2593.3628327722813
Average DIHED: 2935.6501275247506
Average IMPRP: 161.37262722772294
Average ELECT: -5463.052478217821
Average VDW: -1274.3825897029692
###############################

Total run 1010 frames from protein.log.
Average BOND: 952.4267987128711
Average ANGLE: 2593.3628327722813
Average DIHED: 2935.6501275247506
Average IMPRP: 161.37262722772294
Average ELECT: -5463.052478217821
Average VDW: -1274.3825897029692
###############################

Total run 1010 frames from ligand.log.
Average BOND: 952.4267987128711
Average ANGLE: 2593.3628327722813
Average DIHED: 2935.6501275247506
Average IMPRP: 161.37262722772294
Average ELECT: -5463.052478217821
Average VDW: -1274.3825897029692
###############################

delta E(internal): -6642.812386237626
delta E(electrostatic) + deltaG(sol): 5463.052478217821
delta E(VDW) + deltaG(sol): 1274.3825897029692
delta G binding: 94.62268168316405
PS E:\ac\004-MMGBSA>

=====Input for MM calculations with NAMD=====
# ----------input-----

coordinates ionized.pdb
structure ionized.psf
paratypecharmm on
parameters ./par_all36m_prot.prm
parameters ./par_all36_na.prm
parameters ./toppar_water_ions.str
parameters ./par_all35_ethers.prm
parameters ./par_all36_carb.prm
parameters ./par_all36_lipid.prm
parameters ./par_all36_cgenff.prm
parameters ./par_hbond.inp
parameters ./lig.prm
parameters ./ligand.str
#bincoordinates minimized.restart.coor
#binvelocities minimized.restart.vel

# ----------output------

set output output

outputname $output
dcdfile ${output}.dcd
xstFile ${output}.xst
dcdfreq 1000
xstFreq 1000
binaryoutput yes
binaryrestart yes
outputEnergies 1000
restartfreq 1000

# ---------Basic dynamics-------
exclude scaled1-4
1-4scaling 1
COMmotion no
dielectric 1.0

# --------Simulation space partitioning----
switching on
switchdist 9
cutoff 10
pairlistdist 11
MARGIN 3

# --------Multiple time stepping----
firsttimestep 0
timestep 1
stepspercycle 1
rigidbonds all
rigidTolerance 1e-08

# -------Temperature control----
set temperature 300
temperature $temperature; # initial temperature

# -------Langevin Dynamics------
langevin on; # do langevin dynamics
langevinDamping 1; # damping coefficient (gamma) of 1/ps
langevinTemp $temperature; # bath temperature

# ===============================================

PME on
PMEGridSizeX 100
PMEGridSizeY 100
PMEGridSizeZ 100

# doesnt work with just pme

useGroupPressure yes

# with grouppressure, works better, holes still there

LangevinPiston on
LangevinPistonTarget 1.02
LangevinPistonPeriod 150
LangevinPistonDecay 90
LangevinPistonTemp $temperature

# with langevin piston, works!!! But slower.

# ===============================================

# Periodic Boundary conditions
cellBasisVector1 63.62900114059448 0 0
cellBasisVector2 0 71.53299903869629 0
cellBasisVector3 0 0 69.11299800872803
cellOrigin 25.233444213867188 7.527408123016357 21.83285140991211
wrapWater on ;# wrap water to central cell
wrapAll on ;# wrap other molecules too
wrapNearest off

# ---------Scripting
minimize 10000 ;# lower potential energy for 10000 steps
reinitvels $temperature ;# since minimization zeros velocities
run 6000000;

# END
=============Input for GBSA calculations after MM in NAMD=======
#############################################################
## JOB DESCRIPTION ##
#############################################################

# Eq. of KcsA
# embedded in POPC membrane, ions and water.
# Protein released. PME, Constant Pressure/Area.

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################

structure complex.psf
coordinates complex.pdb
outputName complex

set temperature 310

# Continuing a job from the restart files
if {0} {
set inputname LRRK2-VB12-out1
binCoordinates $inputname.restart.coor
binVelocities $inputname.restart.vel ;# remove the "temperature" entry if you use this!
extendedSystem $inputname.restart.xsc
}

firsttimestep 0

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters par_all36_prot.prm
parameters par_all36_cgenff.prm
parameters ligand.str
parameters toppar_water_ions.str

# NOTE: Do not set the initial velocity temperature if you
# have also specified a .vel restart file!
 temperature $temperature

# Periodic Boundary Conditions
# NOTE: Do not set the periodic cell basis if you have also
# specified an .xsc restart file!
if {0} {
cellBasisVector1 98. 0. 0.
cellBasisVector2 0. 98. 0.
cellBasisVector3 0. 0. 96.
cellOrigin -0.0390621498227 -0.0503903478384 0.05063835904
}
# wrapWater on
wrapAll on

### gbsa
GBIS on
solventDielectric 78.5
ionConcentration 0.3
alphaCutoff 15

sasa on
surfaceTension 0.0072

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 16.
switching on
switchdist 15.
pairlistdist 18

# Integrator Parameters
timestep 1 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
nonbondedFreq 2
fullElectFrequency 4
stepspercycle 20

#PME (for full-system periodic electrostatics)
if {0} {
PME yes
PMEGridSizeX 81
PMEGridSizeY 81
PMEGridSizeZ 90
}

# Constant Temperature Control
langevin on ;# do langevin dynamics
langevinDamping 1 ;# damping coefficient (gamma) of 5/ps
langevinTemp $temperature

# Constant Pressure Control (variable volume)
if {0} {
useGroupPressure yes ;# needed for 2fs steps
useFlexibleCell no ;# no for water box, yes for membrane
useConstantArea no ;# no for water box, yes for membrane

langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
langevinPistonPeriod 200.
langevinPistonDecay 50.
langevinPistonTemp $temperature
}

#restartfreq 10000 ;# 1000steps = every 2ps
#dcdfreq 10000
#xstFreq 10000
#outputEnergies 10000
#outputPressure 10000

# Fixed Atoms Constraint (set PDB beta-column to 1)
if {0} {
fixedAtoms on
fixedAtomsFile nottails.fix.pdb
fixedAtomsCol B
fixedAtomsForces on
}

#############################################################
## EXTRA PARAMETERS ##
#############################################################

# Put here any custom parameters that are specific to
# this job (e.g., SMD, TclForces, etc...)

#eFieldOn yes
#eField 0 0 -0.155

#############################################################
## EXECUTION SCRIPT ##
#############################################################

set ts 0

coorfile open dcd complex.dcd

while {[coorfile read] != -1 } {
    firstTimestep $ts
    run 0
    incr ts 1
}
coorfile close

Best regards,
Dr.Shabbir

This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:11 CST