xplor.parseArguments()

import os
print ('process id:'), os.getpid()

numberOfStructures = 100 #number of structure to calculate

import protocol
protocol.initRandomSeed(942)
protocol.initParams("protein")

command = xplor.command
command (""" set message on echo on end """)

import simulation
sim = xplor.simulation

ensembleSize = 2 # 2, 4, 6, 8, 10 members

weights=[1.0/ensembleSize for i in range(ensembleSize)]

pdbTemplate="SCRIPT_%d_STRUCTURE_MEMBER.pdb" % ensembleSize

# create Simulation with ensembleSize ensemble members
from ensembleSimulation import EnsembleSimulation
esim = EnsembleSimulation("ensemble",ensembleSize)
print("ensemble size", esim.size(),"running on",esim.numThreads(),"processors.")

protocol.loadPDB("../initial_structure/ssph1_158-700.pdb",deleteUnknownAtoms=True)

LRR  = "resid 161:393 and segid A"
E3   = "resid 405:622 and segid A"
E2   = "resid 623:697 and segid A"

from potList import PotList
potList = PotList()
crossTerms=PotList('cross terms')
 
from simulationTools import StaticRamp, MultRamp
from simulationTools import InitialParams, FinalParams
rampedParams=[]
highTempParams=[]

from avePot import AvePot
from xplorPot import XplorPot

#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
from repelPotTools import create_RepelPot,initRepel
repel = AvePot(create_RepelPot('repel',
                              selection=AtomSel('not PSEUDO',esim.member())))
potList.append(repel)
rampedParams.append( StaticRamp("initRepel(repel,use14=False)") )
rampedParams.append( MultRamp(.004,4,  "repel.setScale( VALUE)") )

# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""initRepel(repel,
                                               use14=True,
                                               scale=0.004,
                                               repel=1.2,
                                               moveTol=45,
                                               interactingAtoms='name CA'
                                               )""") )

## Selected 1-4 interactions.
from torsionDBPotTools import create_Terminal14Pot
repel14 = AvePot(create_Terminal14Pot('repel14',
                                      selection=AtomSel('not PSEUDO',
                                                        esim.member())))
potList.append(repel14)
highTempParams.append(StaticRamp("repel14.setScale(0)"))
rampedParams.append(MultRamp(0.004, 4, "repel14.setScale(VALUE)"))

##New torsion angle database potential
from torsionDBPotTools import create_TorsionDBPot
tDB = AvePot(create_TorsionDBPot('tDB', system='protein',
                                 selection=AtomSel('known',esim.member())))
potList.append(tDB)
rampedParams.append( MultRamp(.001,2,"tDB.setScale(VALUE)") )

from ensWeightsTools import create_EnsWeights
aWeights=create_EnsWeights("aWeights")

#This energy term implements regularization- so weights don't fluctuate wildly.
potList.append(aWeights)

#set the current ensemble weights
aWeights.setWeights( weights )

#set target weights, used in the regularization energy - by default 1/Ne
aWeights.setTargetWeights( weights )

#set the minimal weight is to 0.05 * 1/Ne  (Ne = enesemble size)
aWeights.setMinFrac(0.025)
#aWeights.setFreedom('fix') # if uncommented, weights will not vary.

#these regularization force constants may need to be adjusted.
rampedParams.append(MultRamp(2500*ensembleSize,20*ensembleSize,
                             "aWeights.setScale(VALUE)"))   
#original scale 50*ensembleSize (force constant at the begining), 0.1*ensembleSize (force constant at the end)


from solnXRayPotTools import create_solnXRayPot
import solnXRayPotTools
xray=create_solnXRayPot('xray',experiment='../restraints/crk2_saxs.dat',
                        numPoints=50, # originally used 30 points
                        maxQ=0.3,
                        normalizeIndex=-3,preweighted=False)

xrayCorrect=create_solnXRayPot('xray-c',experiment='../restraints/crk2_saxs.dat',
                               numPoints=50,
                               maxQ=0.3,
                               normalizeIndex=-3,preweighted=False)

solnXRayPotTools.useGlobs(xray)

xray.setNumAngles(50)
xrayCorrect.setNumAngles(500)
xray.setScale(400)
xray.setCmpType("plain")
potList.append(xray)
crossTerms.append(xrayCorrect)

xray.addEnsWeights( aWeights )
xrayCorrect.addEnsWeights( aWeights )

print (xray.calcEnergy())

from solnScatPotTools import fitParams

#corrects I(q) to higher accuracy, and include solvent contribution corrections
rampedParams.append( StaticRamp("fitParams(xrayCorrect)") )
rampedParams.append( StaticRamp("xray.calcGlobCorrect(xrayCorrect.calcd())") )

#Xplor potentials
bond = AvePot(XplorPot,"BOND")
potList.append(bond)

angl = AvePot(XplorPot,"ANGL")
potList.append(angl)
rampedParams.append( MultRamp(0.4,1,"potList['ANGL'].setScale(VALUE)") )

impr = AvePot(XplorPot,"IMPR")
potList.append(impr)
rampedParams.append( MultRamp(0.1,1,"potList['IMPR'].setScale(VALUE)") )

protocol.massSetup()

from ivm import IVM
dynRigid = IVM()

dynRigid.fix("(%s)" % E3) 
dynRigid.group("(%s)" % LRR) 
dynRigid.group("(%s)" % E2) 

protocol.torsionTopology(dynRigid)

from simulationTools import AnnealIVM
init_t  = 3000 
cool = AnnealIVM(initTemp =init_t,
                 finalTemp=25,
                 tempStep =25, # it was 25
                 ivm=dynRigid,
                 rampedParams = rampedParams)

startPos  = xplor.simulation.atomPosArr()

def calcOneStructure(loopInfo):
    """ this function calculates a single structure, performs analysis on the
    structure, and then writes out a pdb file, with remarks.
    """
    
    xplor.simulation.setAtomPosArr(startPos)
    InitialParams( rampedParams )
  
    # generate a new structure with randomized torsion angles
    from monteCarlo import randomizeTorsions
    randomizeTorsions(dynRigid)
    
    #initial minimization
    protocol.initMinimize(dynRigid,
    potList=potList,
    numSteps=1000,
    printInterval=200)
    dynRigid.run()

    #high-temp dynamics
    protocol.initDynamics(dynRigid,
                          potList=potList,
                          bathTemp=init_t,
                          initVelocities=1,
                          finalTime=800,   # it was 800 
                          numSteps=8000,   # it was 8000 
                          printInterval=200)

    dynRigid.setETolerance( init_t/100 )  #used to det. stepsize. default: t/1000 
    dynRigid.run()
    
    # initialize parameters for cooling loop
    InitialParams( rampedParams )

    # initialize integrator for simulated annealing
    protocol.initDynamics(dynRigid,
                          potList=potList,
                          numSteps=200,       #it was 100 steps or
                          finalTime=.5 ,       # it was 0.2
                          printInterval=100)

    # perform simulated annealing
    cool.run()

    # final torsion angle minimization
    protocol.initMinimize(dynRigid)
    dynRigid.run()

    #do analysis and write structure
#    loopInfo.writeStructure(potList)
    pass

from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
              doWriteStructures=True,
              pdbTemplate=pdbTemplate,
              structLoopAction=calcOneStructure,
              genViolationStats=True,
              averageRegularize=False,
              averagePotList=potList,
              averageCrossTerms=crossTerms,
              averageSortPots=[xray],
              averageTopFraction=0.1, #report on 10% of structs
              averageFitSel="name CA and (%s) " %E3,
              averageContext=FinalParams(rampedParams)).run()