Re: role of gyromagnetic ratio signs in RDC refinement

_______________________________________________
Xplor-nih mailing list
Xplor-nih <at> nmr.cit.nih.gov
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih

Problem running XPLOR

_______________________________________________
Xplor-nih mailing list
Xplor-nih <at> nmr.cit.nih.gov
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih

(no subject)

Hello Charles
After a few attempts to perform dynamics with rigid bodies under RDC
constraints for my multidomain protein, I have to ask again for your help.
I tried to adapt the dock.py script from the eginiput/diffTens directory.
However, instead of getting 10 different structures, all my structures are
exactly the same. IS there something obvious for you that I did wrong?

Below is my script. I hope you will find some time to give a look at it.
Thanks anyway

Hélène Déméné

xplor.requireVersion("2.18")

#
# slow cooling protocol in torsion angle space for protein G. Uses
# NOE, RDC, J-coupling restraints.
#
# this version refines from a reasonable model structure.
#
# CDS 2005/05/10
#

(opts,args) = xplor.parseArguments(["quick"]) # check for command-line typos

quick=False
for opt in opts:
    if opt[0]=="quick":  #specify -quick to just test that the script runs
        quick=True
        pass
    pass

outFilename = "SCRIPT_STRUCTURE.sa"
numberOfStructures=20

if quick:
    numberOfStructures=3
    pass

# protocol module has many high-level helper functions.
#
import protocol
protocol.initRandomSeed(34210)   #explicitly set random seed

#
# annealing settings
#

command = xplor.command

protocol.initParams("protein")

# generate PSF data from sequence and initialize the correct parameters.
#
#from psfGen import seqToPSF
#seqToPSF('protG.seq')
#protocol.initStruct("g_new.psf") # - or from file

# generate a random extended structure with correct covalent geometry
#  saves the generated structure in the indicated file for faster startup
#  next time.
#
#protocol.genExtendedStructure("gb1_extended_%d.pdb" %
#                              protocol.initialRandomSeed())

# or read an existing model
#
protocol.loadPDB("start1.pdb")
xplor.simulation.deleteAtoms("not known")

protocol.fixupCovalentGeom(maxIters=100,useVDW=1)

#
# a PotList contains a list of potential terms. This is used to specify which
# terms are active during refinement.
#
from potList import PotList
potList1 = PotList()
potList0= PotList()
potList2 = PotList()

# parameters to ramp up during the simulated annealing protocol
#
from simulationTools import MultRamp, StaticRamp, InitialParams

rampedParams=[]
highTempParams=[]

# compare atomic Cartesian rmsd with a reference structure
#  backbone and heavy atom RMSDs will be printed in the output
#  structure files
#
#from posDiffPotTools import create_PosDiffPot
#refRMSD = create_PosDiffPot("refRMSD","name CA or name C or name N",
#                            pdbFile='2KUI1.pdb',
#                            cmpSel="not name H*")

# orientation Tensor - used with the dipolar coupling term
#  one for each medium
#   For each medium, specify a name, and initial values of Da, Rh.
#
from varTensorTools import create_VarTensor,addAxisAtoms
media={}
#                        medium  Da   rhombicity
for (medium,Da,Rh,resid) in [
                        ('peg12',   26.6, 0.06,991),
                        ('peg34',   24.9, 0.12,992) ,
         		('pha12',   -11.7, 0.35,995),
                       ('pha34',   -7.3, 0.49,997)
]:

    oTensor = create_VarTensor(medium, axis=addAxisAtoms(resid=resid))
    oTensor.setDaMax(200)
    oTensor.setDa(Da)
    oTensor.setRh(Rh)
    media[medium] = oTensor
    pass

# dipolar coupling restraints for protein amide NH.
#
# collect all RDCs in the rdcs PotList
#
# RDC scaling. Three possible contributions.
#   1) gamma_A * gamma_B / r_AB^3 prefactor. So that the same Da can be used
#      for different expts. in the same medium. Sometimes the data is
#      prescaled so that this is not needed. scale_toNH() is used for this.
#      Note that if the expt. data has been prescaled, the values for rdc
rmsd
#      reported in the output will relative to the scaled values- not the
expt.
#      values.
#   2) expt. error scaling. Used here. A scale factor equal to 1/err^2
#      (relative to that for NH) is used.
#   3) sometimes the reciprocal of the Da^2 is used if there is a large
#      spread in Da values. Not used here.
#
from rdcPotTools import create_RDCPot, scale_toNH
rdcs = PotList('rdc')
for (medium,expt,file,                 scale) in \
    [
     ('pha12','NH','pha12_nh.tbl'       ,1),
   ('pha34','NH','pha34_nh.tbl'       ,1),
    ('peg12','NH' ,'peg12_nh.tbl'       ,1),
     ('peg34','NH','peg34_nh.tbl'       ,1),

     ]:
    rdc = create_RDCPot("%s_%s"%(medium,expt),file,media[medium])

    #1) scale prefactor relative to NH
    #   see python/rdcPotTools.py for exact calculation
    # scale_toNH(rdc) - not needed for these datasets -
    #                        but non-NH reported rmsd values will be wrong.

    #3) Da rescaling factor (separate multiplicative factor)
    # scale *= ( 1. / rdc.oTensor.Da(0) )**2
    rdc.setScale(scale)
    rdc.setShowAllRestraints(1) #all restraints are printed during analysis
    rdc.setThreshold(1.5)       # in Hz
#    rdc.setPotType('square')
    rdcs.append(rdc)
    pass
potList0.append(rdcs)
potList1.append(rdcs)
potList2.append(rdcs)
rampedParams.append( MultRamp(0.05,5.0, "rdcs.setScale( VALUE )") )

# calc. initial tensor orientation
# and setup tensor calculation during simulated annealing
#
from varTensorTools import calcTensorOrientation, calcTensor
for medium in media.keys():
    calcTensorOrientation(media[medium])
    rampedParams.append( StaticRamp("calcTensorOrientation(media['%s'])" %
medium) )
    pass

from simulationTools import  analyze
print analyze(potList0)

# gyration volume term
#
# gyration volume term
#

#Rama torsion angle database
#
from xplorPot import XplorPot
protocol.initRamaDatabase()
potList2.append( XplorPot('RAMA') )
rampedParams.append( MultRamp(.002,1,"potList2['RAMA'].setScale(VALUE)") )

#
# setup parameters for atom-atom repulsive term. (van der Waals-like term)
#
potList2.append( XplorPot('VDW') )
potList1.append( XplorPot('VDW') )
rampedParams.append( StaticRamp("protocol.initNBond()") )
rampedParams.append( MultRamp(0.9,0.8,
                              "command('param nbonds repel VALUE end
end')") )
rampedParams.append( MultRamp(.004,4,
                              "command('param nbonds rcon VALUE end end')") )
# nonbonded interaction only between CA atoms
highTempParams.append( StaticRamp("""protocol.initNBond(cutnb=100,
                                                        rcon=0.004,
                                                        tolerance=45,
                                                        repel=1.2,
                                                        onlyCA=1)""") )

potList2.append( XplorPot("BOND") )
potList2.append( XplorPot("ANGL") )
potList2['ANGL'].setThreshold( 5 )
rampedParams.append( MultRamp(0.4,1,"potList2['ANGL'].setScale(VALUE)") )
potList2.append( XplorPot("IMPR") )
potList2['IMPR'].setThreshold( 5 )
rampedParams.append( MultRamp(0.1,1,"potList2['IMPR'].setScale(VALUE)") )

# Give atoms uniform weights, except for the anisotropy axis
#

# IVM setup
#   the IVM is used for performing dynamics and minimization in torsion-angle
#   space, and in Cartesian space.
#

# initially minimize in Cartesian space with only the covalent constraints.
#   Note that bonds, angles and many impropers can't change with the
#   internal torsion-angle dynamics
#   breaks bonds topologically - doesn't change force field
#
#dyn.potList().add( XplorPot("BOND") )
#dyn.potList().add( XplorPot("ANGL") )
#dyn.potList().add( XplorPot("IMPR") )
#
#dyn.breakAllBondsIn("not resname ANI")
#import varTensorTools
#for m in media.values():
#    m.setFreedom("fixDa, fixRh")                 #fix tensor parameters
#    varTensorTools.topologySetup(dyn,m) #setup tensor topology

#protocol.initMinimize(dyn,numSteps=1000)
#dyn.run()

# reset ivm topology for torsion-angle dynamics
#

from ivm import IVM
from ivm import PublicIVM
dyn_fix = IVM()
dyn_fix.reset()
dyn_fix.group( 'resid 358:421' )
dyn_fix.group( 'resid 426:488' )
dyn_fix.group( 'resid 496:552' )
dyn_fix.group( 'resid 556:625' )
import varTensorTools

for m in media.values():
    m.setFreedom("fixDa, fixRh")        #fix tensor Rh, Da, vary orientation
#    m.setFreedom("varyDa, varyRh")      #vary tensor Rh, Da, vary
orientation
    pass

varTensorTools.topologySetup(dyn_fix,list=media.values())

protocol.torsionTopology(dyn_fix,oTensors=media.values())

dyn_free_sch = IVM()
dyn_free_sch.reset()
dyn_free_sch.group( 'resid 358:421 and \
                             (name CA or  name C or name N or name O or
name HN)' )
dyn_free_sch.group( 'resid 426:488 and \
                             (name CA or  name C or name N or name O or
name HN)')

dyn_free_sch.group( 'resid 496:552 and  \
(name CA or  name C or name N or name O or name HN)')
dyn_free_sch.group( 'resid 556:625 and \
(name CA or  name C or name N or name O or name HN)')

for m in media.values():
    m.setFreedom("fixDa, fixRh")        #fix tensor Rh, Da, vary orientation
#    m.setFreedom("varyDa, varyRh")      #vary tensor Rh, Da, vary
orientation
    pass

varTensorTools.topologySetup(dyn_free_sch,list=media.values())

protocol.torsionTopology(dyn_free_sch,oTensors=media.values())

# object which performs simulated annealing
#
from simulationTools import AnnealIVM
init_t  = 3000.     # Need high temp and slow annealing to converge
cool = AnnealIVM(initTemp =init_t,
                 finalTemp=10,
                 tempStep =12.5,
                 ivm=dyn_fix,
                 rampedParams = rampedParams)

# Give atoms uniform weights, except for the anisotropy axis
#
from atomAction import SetProperty
import varTensorTools
AtomSel("not resname ANI").apply( SetProperty("mass",100.) )
AtomSel("all            ").apply( SetProperty("fric",10.) )

def accept(potList):
    """
    return True if current structure meets acceptance criteria
    """
    if potList['IMPR'].violations()>1:
        return False

    return True

def calcOneStructure(loopInfo):
    """ this function calculates a single structure, performs analysis on the
    structure, and then writes out a pdb file, with remarks.
    """

    from simulationTools import  analyze
    print analyze(potList0)

    # initialize parameters for high temp dynamics.

    InitialParams( rampedParams )
    # high-temp dynamics setup - only need to specify parameters which
    #   differfrom initial values in rampedParams
    InitialParams( highTempParams )

    # initial minimization

    protocol.initMinimize(dyn_fix, potList=potList0,
                            printInterval=50)

    dyn_fix.run()

    dyn_fix.run()

    protocol.initMinimize(dyn_fix, potList=potList0,
                            printInterval=50)

    dyn_fix.run()

    dyn_fix.run()

    protocol.initMinimize(dyn_fix, potList=potList0,
                            printInterval=50)

    dyn_fix.run()

    dyn_fix.run()

    dyn_fix.run()

#    k2=0
#    while  k2 < 0:
    # initialize parameters for cooling loop
    InitialParams( rampedParams )

     # initialize integrator for simulated annealing
     #
    protocol.initDynamics(dyn_fix,
                          potList=potList2,
                          numSteps=1000,       #at each temp: 100 steps or
                          finalTime=1 ,       # 1 ps, whichever is less
                          printInterval=100)

      # perform simulated annealing
      #
    cool.run()
#    k2=k2+1
#    print k2

    # final torsion angle minimization
    #
    protocol.initMinimize(dyn_free_sch,potList=potList2,
                          printInterval=50)
    dyn_fix.run()

    # final all- atom minimization
    #

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

from simulationTools import StructureLoop, FinalParams
StructureLoop(numStructures=numberOfStructures,
              pdbTemplate=outFilename,
              structLoopAction=calcOneStructure,
              genViolationStats=1,
              averagePotList=potList2,
              averageSortPots=[potList2['BOND'],potList2['ANGL'],potList2['IMPR'],
                               rdcs],
#              averageCrossTerms=refRMSD,
              averageTopFraction=0.5, #report only on best 50% of structs
              averageAccept=accept,   #only use structures which pass
accept()
              averageContext=FinalParams(rampedParams),
              averageFilename="SCRIPT_ave.pdb",    #generate regularized
ave structure
              averageFitSel="name CA",
              averageCompSel="not resname ANI and not name H*"     ).run()

--

-- 
Passerelle antivirus CBS

Re: (no subject)

Hello Hélène--

> After a few attempts to perform dynamics with rigid bodies under RDC
> constraints for my multidomain protein, I have to ask again for your help.
> I tried to adapt the dock.py script from the eginiput/diffTens directory.
> However, instead of getting 10 different structures, all my structures are
> exactly the same. IS there something obvious for you that I did wrong?
> 

Do did leave out the position/orientation randomization step. This
should be added to the beginning of calcOneStructure:

    from atomAction import randomizeDomainPos
    randomizeDomainPos( 'resid 1:249', deltaPos=45 )

where the atom selection specified one of the domains.

best regards--
Charles

--
Charles Schwieters     email:   Charles <at> Schwieters.org
                       www:     http://schwieters.org/cds
phone: (301) 402-4914  PGP key: http://schwieters.org/cds/pgp.txt

_______________________________________________
Xplor-nih mailing list
Xplor-nih <at> nmr.cit.nih.gov
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih

using chromophore residues from toph19.chromo to generate psf

Dear all,

I would like to calculate a structure containing a bacteriochlorophyll  
a (bchla) ligand. I have found what I believe is a defintion of  
parameters for this ligand in the file "toph19.chromo". In this file I  
would think that the ligand parameters are wrapped up in the
"RESIdue BCL" and following statements

As the first step I tried to generate the psf file from sequence  
templates using the BCL residue for the ligand. A problem occurs when  
I try loading the parameters, maybe I do in the wrong place or  
context, hope your can help me, here is my test script:

parameter  <at> TOPPAR:protein.par  <at> TOPPAR:toph19.chromo end

segment
    name="A   "
    SETUP=TRUE
    chain
        <at> TOPPAR:toph11.pep
sequence
HIS BCL
end
end
end

write psf output=test.psf end

This leads to parsing errors regardsless of if I put the loading of  
BCL parameters within the chain block. The parsing errors seams to be  
related to the AUTOGENERATE keyword, which is not recognized, in the  
toph19.chromo. Here is that top part of the content of toph19.chromo:

remarks   toph19.chromo (CHROMOPHORE RESIDUE'S )
remarks   ================================================
remarks   Authors: Herbert Treutlein and Axel Brunger
remarks   warning: this file is under development

set echo=false end

AUTOGENERATE  ANGLES=TRUE  END      {* but see the two omitted angles below *}
{*============================         ==================================== *}

{* protein default masses *}
MASS   H      1.00800! hydrogen which can h-bond to neutral atom
MASS   HC     1.00800!    ="=     ="=     ="=    to charged atom

Hope you can help to point out where I am wrong.

Furthermore, the bchla in our structure is bound to a histidine. Would  
it be advisable to use the PBCL patch? That means that I use PBCL in  
place of HIS and BCL?

best regards,

Jakob

I am using xplor-nih 2.26

--

-- 
Jakob Toudahl Nielsen, post doc
Laboratory for Biomolecular NMR Spectroscopy
inSPIN, Center for Insoluble Protein Structures
Department of Chemistry, University of Aarhus
Langelandsgade 140, DK-8000 Aarhus C, Denmark
Phone: +458942 5528 or +452993 8502 (cell)

Re: using chromophore residues from toph19.chromo to generate psf

Dear all,

sory, I realize that I need the "param19.chromo" as the parameter  
file. However, this time I still get errors, it seams that the BCL  
residue is not recognized. Here is output from the psf generation,  
which fails to generate parameters for BCL.

...
  CHAIN>sequence
  SEQUENCE>HIS BCL
  MAPIC: Atom numbers being modified
  %GENRES-ERR: residue BCL  not found in topology file
  %PATCH-ERR: bond -C    +N    not found
  %PATCH-ERR: angle -CA   -C    +N    not found
  %PATCH-ERR: angle -O    -C    +N    not found
... more PATCH-ERRs...

hope that someone can help to find out what I do wrong.

best regards,

Jakob

Citing Jakob Toudahl Nielsen <jtn <at> chem.au.dk>:

> Dear all,
>
> I would like to calculate a structure containing a bacteriochlorophyll
> a (bchla) ligand. I have found what I believe is a defintion of
> parameters for this ligand in the file "toph19.chromo". In this file I
> would think that the ligand parameters are wrapped up in the
> "RESIdue BCL" and following statements
>
> As the first step I tried to generate the psf file from sequence
> templates using the BCL residue for the ligand. A problem occurs when
> I try loading the parameters, maybe I do in the wrong place or
> context, hope your can help me, here is my test script:
>
> parameter  <at> TOPPAR:protein.par  <at> TOPPAR:toph19.chromo end
>
> segment
>     name="A   "
>     SETUP=TRUE
>     chain
>         <at> TOPPAR:toph11.pep
> sequence
> HIS BCL
> end
> end
> end
>
> write psf output=test.psf end
>
> This leads to parsing errors regardsless of if I put the loading of
> BCL parameters within the chain block. The parsing errors seams to be
> related to the AUTOGENERATE keyword, which is not recognized, in the
> toph19.chromo. Here is that top part of the content of toph19.chromo:
>
> remarks   toph19.chromo (CHROMOPHORE RESIDUE'S )
> remarks   ================================================
> remarks   Authors: Herbert Treutlein and Axel Brunger
> remarks   warning: this file is under development
>
> set echo=false end
>
> AUTOGENERATE  ANGLES=TRUE  END      {* but see the two omitted  
> angles below *}
> {*============================          
> ==================================== *}
>
> {* protein default masses *}
> MASS   H      1.00800! hydrogen which can h-bond to neutral atom
> MASS   HC     1.00800!    ="=     ="=     ="=    to charged atom
>
> Hope you can help to point out where I am wrong.
>
>
> Furthermore, the bchla in our structure is bound to a histidine. Would
> it be advisable to use the PBCL patch? That means that I use PBCL in
> place of HIS and BCL?
>
>
> best regards,
>
> Jakob
>
> I am using xplor-nih 2.26
>
> --
> Jakob Toudahl Nielsen, post doc
> Laboratory for Biomolecular NMR Spectroscopy
> inSPIN, Center for Insoluble Protein Structures
> Department of Chemistry, University of Aarhus
> Langelandsgade 140, DK-8000 Aarhus C, Denmark
> Phone: +458942 5528 or +452993 8502 (cell)
> _______________________________________________
> Xplor-nih mailing list
> Xplor-nih <at> nmr.cit.nih.gov
> http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih
>

--

-- 
Jakob Toudahl Nielsen, post doc
Laboratory for Biomolecular NMR Spectroscopy
inSPIN, Center for Insoluble Protein Structures
Department of Chemistry, University of Aarhus
Langelandsgade 140, DK-8000 Aarhus C, Denmark
Phone: +458942 5528 or +452993 8502 (cell)

Re: Xplor-nih Digest, Vol 103, Issue 4

_______________________________________________
Xplor-nih mailing list
Xplor-nih <at> nmr.cit.nih.gov
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih

Creating a psf from a pdb using psfGen

Hi All,

I am attempting to use a new python script to perform some PRE 
calculations and I am using a structure containing 6 spin labels with 
four orientations attached to a haem containing protein, and when I run 
psfGen to obtain the psf file I get the following traceback:

/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py:220: 
DeprecationWarning: raising a string exception is deprecated
  raise "residue %s not of type %s in %s" % (res,ret,`seq`)
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File 
"/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/bin/pdb2psf", 
line 47, in <module>
    pdbToPSF(open(infile).read())
  File 
"/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py", 
line 588, in pdbToPSF
    deprotonateHIS=False,customRename=customRename)
  File 
"/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py", 
line 355, in seqToPSF
    if seqType=='auto': seqType = deduceSeqType(seq)
  File 
"/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py", 
line 220, in deduceSeqType
    raise "residue %s not of type %s in %s" % (res,ret,`seq`)
 Is there a way to do this in pyXplor or do I need to use the 
old-fashioned xplor to achieve this?

Alternatively, is there a way to re-number residues in a psf to 
correspond with those in a pdb file?

Regards

Simon

--

-- 
Simon P Skinner
Protein Chemistry Group
Leiden Institute of Chemistry, Universiteit Leiden
Phone: +31 71 527 6089 / Fax: +31 71 527 4349
E-mail : skinnersp <at> chem.leidenuniv.nl 

Re: Creating a psf from a pdb using psfGen

Hello Simon--

> I am attempting to use a new python script to perform some PRE
> calculations and I am using a structure containing 6 spin labels with
> four orientations attached to a haem containing protein, and when I run
> psfGen to obtain the psf file I get the following traceback:
>
> /afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py:220:
> DeprecationWarning: raising a string exception is deprecated
>   raise "residue %s not of type %s in %s" % (res,ret,`seq`)
> Traceback (most recent call last):
>   File "<string>", line 1, in <module>
>   File
> "/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/bin/pdb2psf",
> line 47, in <module>
>     pdbToPSF(open(infile).read())
>   File
> "/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py",
> line 588, in pdbToPSF
>     deprotonateHIS=False,customRename=customRename)
>   File
> "/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py",
> line 355, in seqToPSF
>     if seqType=='auto': seqType = deduceSeqType(seq)
>   File
> "/afs/gorlaeus.net/lic/metprot/software/xplor-nih/2.22/python/psfGen.py",
> line 220, in deduceSeqType
>     raise "residue %s not of type %s in %s" % (res,ret,`seq`)
>  Is there a way to do this in pyXplor or do I need to use the
> old-fashioned xplor to achieve this?

You might check out eginput/pre/generate.py. Here's another example
which uses a combination of scripts. This will be distributed in the
next release of Xplor-NIH

  http://nmr.cit.nih.gov/pub/pre2.tar.gz

hope it helps--
Charles

Charles Schwieters     email:   Charles <at> Schwieters.org
                       www:     http://schwieters.org/cds
phone: (301) 402-4914  PGP key: http://schwieters.org/cds/pgp.txt

_______________________________________________
Xplor-nih mailing list
Xplor-nih <at> nmr.cit.nih.gov
http://dcb.cit.nih.gov/mailman/listinfo/xplor-nih