Dear all,
I'm now refining a G-quadruplex structure using a modified refine.py
based on the script in dir egiput/dna_refi of xplor-NIH. The input data
is only noe constraints.
My question is if I use Lennard-Jones and electrostatics and add the lines
in refine.py:
rampedParams.append( StaticRamp("""xplor.command('''param nbonds
atom
repel=0.8
wmin=0.01
nbxmod=3
cutnb=11.5
cctonnb=9.5
ctofnb=10.5
tolerance=0.5
rdie
vswitch
switch
end end''')"""))
the energies for noe,vdw, bond,angle etc. are very high and there're lots
of noe violations(> 0.5) which are zero in the initial calculation. But if
use :
protocol.initNBond(cutnb=4.5)
potList.append( XplorPot("VDW") )
potList.append( XplorPot("elec") )
rampedParams.append( MultRamp(0.9,0.78,
"xplor.command('param nbonds repel VALUE end
end')") )
rampedParams.append( MultRamp(.004,4,
"xplor.command('param nbonds rcon VALUE end
end')") )
The output is fine.But the program will give the message:
InternalDynamics::step: large timestep detected. Halving.
Would you please give any explanations and suggestions?
Thanks in advance!
Best
Dong
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
below is part of refine.py:
protocol.initNBond(cutnb=4.5)
potList.append( XplorPot("VDW") )
potList.append( XplorPot("elec") )
rampedParams.append( MultRamp(0.9,0.78,
"xplor.command('param nbonds repel VALUE end
end')") )
rampedParams.append( MultRamp(.004,4,
"xplor.command('param nbonds rcon VALUE end
end')") )
#highTempParams.append( StaticRamp("""protocol.initNBond(cutnb=100,
# tolerance=45,
# #repel=1.2,
# #selStr='name p'
# )"""))
for name in ("bond","angl","impr"):
potList.append( XplorPot(name) )
pass
rampedParams.append( MultRamp(0.4,1.0,"potList['ANGL'].setScale(VALUE)"))
rampedParams.append( MultRamp(0.1,1.0,"potList['IMPR'].setScale(VALUE)"))
from ivm import IVM
dyn = IVM()
protocol.initDynamics(dyn,potList=potList)
protocol.torsionTopology(dyn)
# Give atoms uniform weights, except for the anisotropy axis
from atomAction import SetProperty
AtomSel("not resname ANI").apply( SetProperty("mass",100.) )
#varTensorTools.massSetup(media.values(),300)
AtomSel("all ").apply( SetProperty("fric",10.) )
##
## minc used for final cartesian minimization
##
from selectTools import IVM_groupRigidSidechain
minc = IVM()
protocol.initMinimize(minc,potList=potList)
IVM_groupRigidSidechain(minc)
protocol.cartesianTopology(minc,"not resname ANI")
#varTensorTools.topologySetup(minc,media.values())
init_t1 = 20000
init_t2 = 3000
from simulationTools import AnnealIVM
anneal1= AnnealIVM(initTemp =init_t1,
finalTemp=init_t2,
tempStep =500,
ivm=dyn,
rampedParams = rampedParams)
anneal2= AnnealIVM(initTemp =init_t3,
finalTemp=25,
tempStep =25,
ivm=dyn,
rampedParams = rampedParams)
# initialize parameters for initial minimization.
InitialParams( rampedParams )
# initial minimization
InitialParams( highTempParams )
protocol.initMinimize(dyn,
potList=[potList['CDIH'],potList['IMPR']],
numSteps=50)
dyn.run()
####################################################################
protocol.initMinimize(dyn,
potList=potList,
numSteps=1000)
minc.run()
def calcOneStructure(loopInfo):
InitialParams( rampedParams )
InitialParams( highTempParams )
protocol.initDynamics(dyn,
initVelocities=1,
bathTemp=init_t2,
potList=potList,
finalTime=50)
dyn.setETolerance( init_t2/100)
dyn.run()
InitialParams( rampedParams )
protocol.initDynamics(dyn,
finalTime=0.5,
numSteps=0,
#eTol_minimum=0.001
)
anneal1.run()
anneal2.run()
#
# torsion angle minimization
#
protocol.initMinimize(dyn,numSteps=5000)
dyn.run()
##
##all atom minimization
##
protocol.initMinimize(minc,potList=potList,numSteps=3000)
minc.run()
#
# perform analysis and write structure
loopInfo.writeStructure(potList,crossTerms)
pass
from simulationTools import StructureLoop
StructureLoop(numStructures=numberOfStructures,
startStructure=startStructure,
structLoopAction=calcOneStructure,
pdbTemplate=outFilename,
genViolationStats=1,
averageFilename="average_min.pdb",
#averageFitSel="",
averageRefineSteps=15,
averageTopFraction=0.3,
averagePotList=potList).run()
