Dear Charles and xplor-users,

I am trying to use ensemble calculation for unstructured proteins.  However, the output structures are too compact than expected, so I did the following test:

The simulated annealing python input script from the PRE distribution was used, but removed all experimental restraint terms.  Only "angl", "vdw", "bond", "impr", and "rama" were  left.

And I added or modified the following lines in proper positions of this script for ensemble calculation:
########################################
from ensembleSimulation import EnsembleSimulation

ensembleSize=1
esim = EnsembleSimulation("ensemble",ensembleSize)

hitemp_potList = PotList()
xpTerms = ["BOND", "ANGL", "IMPR", "VDW", "RAMA"]
for xn in  xpTerms :
   hitemp_potList.add( AvePot(XplorPot, xn)  )
   pass

dyn  = IVM(esim)
########################################
If the ensemble size equals to 1,  the results are quite reasonable:  I calculated 500 times (numStructures=500) and analyzed their average radius of gyration, their sizes are widely distributed with a averaged value close to  expected random coil size(in my case, 76 residues with radius of gyration around 40 A).

However, if the ensemble size is set to 7, the output structures are much compact, the radius of gyration is only around 26A.  Something wrong with the ensemble calling?  Thanks in advanced!!

Best wishes,

Jie-rong

Here is the whole script I used for ensemble calculation.
########################
#
#
command = xplor.command
from jCoupPot import JCoupPot
from noePot import NOEPot
from xplor import select
from xplorPot import XplorPot
from rdcPotTools import *
from pdbTool import *
from atomAction import *
from selectTools import *
from simulationTools import *
from ivm import IVM
import string
import protocol

from ensembleSimulation import EnsembleSimulation

simWorld.setRandomSeed( 5764 )

command("""
 structure <at> ubq1conf.psf end
 param <at> ./TopPar/parallhdg_new.pro
       <at> ./TopPar/parnah1er1_mod_new.inp
       <at> ./TopPar/par_axis_3.pro
       <at> ./TopPar/edtaThymine.par
       <at> ./TopPar/manganese.par
 end
 set echo = on message = on end
""")
protocol.initParams(('./TopPar/cysp.par'))
command("""
 coor <at> random1conf1.pdb
 coor copy end
""")

ensembleSize=7
esim = EnsembleSimulation("ensemble",ensembleSize)

tmps = open("./TopPar/proShortRange.list").read()
proShortRange = tmps.split(None)

#
#Rama torsion angle databasse
#

protocol.initRamaDatabase()

#
# setup parameters for nbond term.
#

command("""
   parameters  ! just repulsive term
    nbonds
      atom
      nbxmod 3  wmin 0.01 cutnb 4.5 tolerance 0.5
      rexp 2 irex 2 ctonnb=2.99 ctofnb=3.
      repel= 0.5 rcon= 1.0 ! changed later
    end
   end
""")

#
# annealing setting
#

ini_rad  = 0.9  ; fin_rad = 0.78
ini_con = 0.004 ; fin_con = 4.0
ini_ang = 0.4   ; fin_ang = 1.0
ini_imp = 0.1   ; fin_imp = 1.0

ini_rama =  0.002 ; fin_rama =  1.0

hot_cdi =  10.0   ; cool_cdi = 200.0

hitemp_potList = PotList()
xpTerms = ["BOND", "ANGL", "IMPR", "VDW", "RAMA"]
for xn in  xpTerms :
   hitemp_potList.add( AvePot(XplorPot, xn)  )
   pass

#set mass and fric
AtomSel("not (resname ANI or resname TAU)").apply(SetProperty("mass", 100.0))
AtomSel("resname ANI or resname TAU ").apply(SetProperty("mass", 300.0))
AtomSel("all").apply(SetProperty("fric", 10.0))

#
# IVM setup
#

dyn  = IVM(esim)
dyn.reset() ; dyn.potList().removeAll()

#break ribose rings
IVM_breakRiboses(dyn)

#break proline rings
IVM_breakProlines(dyn)

IVM_groupRigidSidechain(dyn)
dyn.autoTorsion()

def setConstraints(k_ang,k_imp):
    nconf = 1
    scaleVdw = 1.0 / float( nconf )
    command("""
      constraints
         inter = (segid "C1") (segid "C1")
         weights * 1 angl %f impr %f vdw %f end
      end""" % ( k_ang, k_imp, scaleVdw) )
    return

def checkAndCorrect(dyn, pene, ene, tmstp, ni, maxEner):
    # successful or not?
    #
    badIter = 0
    if ene > maxEner :
      badIter = 1
      print "energy is too high (%7.1f)" % ene
    if ene > pene and ene/pene > 10.0 :
      badIter = 1
      print "energy is higher than ten times of previous value (current: %7.1f, previous: %7.1f)" % (ene, pene)
    if tmstp < 1.0e-7 : #oritinal 1.0e-4
      print "timestep is too short (%e ps)" % tmstp
      badIter = 1

    # for successful case
    #  - copy the coordinates to comparison set
    #
    if badIter == 0 :
      command("""
        coor copy end ! copy good coordinates
      """)
      return 1
  
    # for unsuccessful case
    #  - copy back the previous coordinates 
    #  - randomize initial velocities

    command("""
      coor swap end ! call previous good coordinates
      coor copy end ! <- to erase bad one
    """)

    if ni > 3 :
       print "Too many corrections. Exiting programs"
       exit(1)

    bath = dyn.bathTemp()
    randomizeVelocities( bath )
    dyn.resetReuse()  #necessary to get randomized velocities
    return 0

def dynamicsWithCheck(s, timestep, bath, finalTime, maxcycles):
    s.calcEnergy()
    ok = 0
    ni = 0
    while ok == 0 :
         pene   = s.Epotential()
         s.setNumSteps(maxcycles)
         s.setStepsize( timestep )
         s.setBathTemp( bath )
         s.setFinalTime( finalTime )
         s.setPrintInterval( 500 ) #original 10
         s.setETolerance( bath/1000 )
         s.setResetCMInterval(10) 
         s.setResponseTime(20)
         #dyn.setVerbose(0xffff)
         dyn.setMaxDeltaE(10000)
         #dyn.setAdjustStepsize(0)
         s.run()
         ene    = s.Epotential()
         tmstp  = s.stepsize()
         ok = checkAndCorrect(s, pene, ene, tmstp, ni, 100000) 
         ni = ni + 1
         pass
    return

def setProRamaForce(k_rama):
    trace.suspend()
    for cn in (proShortRange):
       command("rama class %s force %f end " % (cn, k_rama))
       pass
    trace.resume()
    print "(setProRamaForce) current k_rama: ", k_rama
    return

def structLoopAction(loopInfo):

    #
    # read initial structure
    #
    command("""
     coor init end
     coor <at> random1conf1.pdb
     coor copy end
    """)

    dyn.init()
    #
    # regularization of covalent geometry
    #                      
    command("""
      flags exclude * include bond angl impr end
      mini powell nstep 1000 nprint 50 end
    """)

    #
    # high temp dynamics with vdw
    #

    # potential terms to use
    dyn.setPotList( hitemp_potList)

    setProRamaForce( ini_rama )
    command("restraints dihed scale %f end" % hot_cdi )

    # VDW setting for high temperature
    command("""
      parameters
        nbonds
          atom
          nbxmod 3
          wmin  =   0.01  cutnb =   100   tolerance 45
          repel= 1.2   rexp   =  2   irex   =  2  rcon=1.0 
        end
      end
      constraints
       inter = (not (name ca or name p or name mn)) (all)
       weights * 1 angl 0.4 impr 0.1 vdw 0 elec 0 end
       inter = (name ca or name p or name mn) (name ca or name p or name mn)
       weights * 1 angl 0.4 impr 0.1 vdw 1.0 end
      end
    """)

    # For PRE-calc in SB mode
 
    # dynamics
    timestep  = 0.002
    bath = 3000.0
    dyn.setStepType("pc6")
    randomizeVelocities( bath )
    dyn.resetReuse()  #necessary to get randomized velocities
    dyn.calcEnergy()
    for i_high in range(0,10) :
       print " high temp dynamics -- subcycle No. %d" % i_high
       dynamicsWithCheck(dyn, timestep, bath, 1.0, 500)
       #print pre.showRestraints(0)
       pass


    #
    # cooling
    #

    #number of cycles
    numSteps = 240

    global k_ang, k_imp

    # VDW setting for cooling
    command("""
     parameters
      nbonds
       atom
       nbxmod 3 
       wmin 0.01 
       cutnb 4.5 
       tolerance 0.5
       rexp 2 irex 2
       ctonnb=2.99 ctofnb=3.
       repel=0.5  ! changed later
       rcon=1.0   ! changed later
      end
     end
    """)
    setConstraints(ini_ang, ini_imp)

    # setting force constants and etc.
    k_rama = MultRamp(ini_rama,fin_rama,"setProRamaForce( VALUE )")
    k_ang  = MultRamp(ini_ang,fin_ang)
    k_imp  = MultRamp(ini_imp,fin_imp)
    radius = MultRamp(ini_rad,fin_rad, "command('param nbonds repel VALUE end end')")
    k_vdw  = MultRamp(ini_con,fin_con, "command('param nbonds rcon VALUE end end')")

    rampedParams = [k_rama,radius,k_vdw,k_ang,k_imp]

    for k in rampedParams :
       k.init(numSteps)
       pass
    command("restraints dihed scale %f end" % cool_cdi )

    #dynamics
    initTemp  = 3000.0
    finalTemp =  25.0  #original 25
    coolTimestep = 0.002
    tempStep  = (initTemp-finalTemp)/float(numSteps)
    finalTime = 1.0      #original 0.2
    #maxCycle =  100     #<==original value
    maxCycle =  500
    bathTemp = initTemp;
    for i_cool in range(0,numSteps):
       bathTemp -= tempStep
       for k in rampedParams :
            k.update()
            pass
       setConstraints(k_ang.value(), k_imp.value())
       print "cooling cycle - No.%d-" % i_cool
       print "current bathTemp = %7.2f" % bathTemp
       dynamicsWithCheck(dyn, coolTimestep, bathTemp, finalTime, maxCycle)
       pass

    #
    # Powell minimization in torsion angle space
    #
    dyn.setStepType("powell")
    dyn.setNumSteps(20000)
    dyn.setMaxCalls(20000)
    dyn.setETolerance(1.0e-7)
    dyn.setGTolerance(0.01)
    dyn.setDEpred(1.0)
    dyn.run()

    #
    # Powell minimization in Cartesian space
    #
    xplor.simulation.potList().removeAll()

    command("""
    flags exclude *
          !include bond angle impr vdw rama orie plan noe cdih coup carb sani xdip coll scripting
      include bond angle impr vdw rama
    end
    mini powell nstep 2000 nprint 10 end
    """)

    dyn.resetReuse()
    dyn.init()
    dyn.calcEnergy()
   
    outPDB_template = "./RandomConfNOPRE/%d_MEMBER_%d.pdb" % (ensembleSize, (loopInfo.count+1))
    outFile = PDBTool( outPDB_template )

    en = {}
    eAll  = dyn.Epotential()

    rms = {}
    violations = {}
    for (term, threshold) in (("BOND",0.05),
                              ("ANGL", 5.0),
                              ("IMPR", 5.0),
                              ("NOE",  0.5),
                              ("CDIH", 5.0)
                  ):
        result = command("print thresh %f %s"% (threshold,term),
                         ("RESULT","VIOLATIONS&q uot;))
        rms[term] = float(result[0])
        violations[term] = int(result[1])
        pass


    outFile.addRemark("OVERALL ENERGY: %7.2f " % eAll)
    bars = "----------------------------------------------------------"
    outFile.addRemark(bars)
    outFile.addRemark("RMS")
    outFile.addRemark("bond: %7.4f, angl: %7.4f, impr: %7.4f" % (rms["BOND"], rms["ANGL"], rms["IMPR"]))

    outFile.write()
    pass

StructureLoop(numStructures=50,startStructure=0 ,
              structLoopAction=structLoopAction).run()

Dear Charles,

Thanks for your suggestion.  I tried to understand this order.py script from gb3_ensemble, but there are still some parts that still confused me.

The biggest difference I can tell between the script I posted in previous e-mail and order.py is that  there is a Powell minimization following the dynamics at each temperature during the cooling steps.  Is this the feature that you suggested "can't mix EnsembleSimulaition with XPLOR minimization"?  I don't understand this.  Could you explain more?  Or some reference suggested to read?

You also mentioned "the scaling of energy results in a time scaling, such that larger ensembles should use longer times for dynamics".  Does it mean that if larger ensemble size is used, the dynamics time should set longer manually (for example, set finalTime 1 ps to 10 ps )?  Or does it mean the "delta_t" between dynamic steps become shorter when the size of ensemble is larger, and it takes longer (more computer time, more steps) to reach finalTime?

In addition, any difference between potList.add, and potList.append, since both are used in order.py?  Thanks a lot!!

Best wishes,

Jie-rong