diff --git a/PyAttract/attract.py b/PyAttract/attract.py index 4d2016a..237e5e2 100755 --- a/PyAttract/attract.py +++ b/PyAttract/attract.py @@ -11,8 +11,6 @@ import bz2 #for compression of Ligand and receptor data import base64 #compressed ligand and receptor as base64 strings - - def surreal(i): return i @@ -33,6 +31,17 @@ def PrintVect(vect): print vect[i], " | ", print '' +def printWeights(w, e): + print "### WEIGHTS BEGIN" + for region in xrange(0, len(w)): + sum = 0 + print "### Start Region %d" %(region+1) + for copy in xrange(0, len(w[region])): + sum += w[region][copy] + print "WEIGHT REGION %d COPY %d = %f (Energy = %f)" %(region+1, copy+1, w[region][copy], e[region][copy]) + print "WEIGHT SUM OF WEIGHTS = %f" %sum + print "### End Region %d" %(region+1) + print "### WEIGHTS END" class Rotation: @@ -85,6 +94,9 @@ def __init__(self): def __iter__(self): return self._rot.__iter__() + def __getitem__(self, i): + return self._rot.__getitem__(i) + class Translation: @@ -185,19 +197,18 @@ def checkFile(name, comment): ## MAIN ATTRACT PROGRAM # ########################### from optparse import OptionParser -parser = OptionParser(usage="%prog -r receptor_file -l ligand_file [-h] [-s] [-t] [--ref]") +parser = OptionParser(usage="%prog -r receptor_file -l ligand_file [-h] [-s] [-trans] [-rot] [--ref]") parser.add_option("-r", "--receptor", action="store", type="string", dest="receptor_name", help="name of the receptor file") parser.add_option("-l", "--ligand", action="store", type="string", dest="ligand_name", help="name of the ligand file") parser.add_option("-s", "--single", action="store_true", dest="single", default=False, help="single minimization mode") parser.add_option("--ref", action="store", type="string", dest="reffile", help="reference ligand for rmsd" ) -parser.add_option("-t", "--translation", action="store", type="int", dest="transnb", help="translation number (distributed mode) starting from 0 for the first one!") -parser.add_option("--start1", action="store_true", default=False, dest="start1", help="(only useful with -t), use 1 for the first translation point") +parser.add_option("--trans", "--translation", action="store", type="int", dest="transnb", help="translation number (distributed mode).") +parser.add_option("--rot", "--rotation", action="store", type="int", dest="rotnb", help="rotation number (distributed mode).") +parser.add_option("--start1", action="store_true", default=False, dest="start1", help="(only useful with -trans), use 1 for the first translation point") +parser.add_option("--mcop", action="store_true", default=False, dest="regions", help="mcop option for multi-copy rigid-body docking") (options, args) = parser.parse_args() -#receptor_name=args[0] -#ligand_name=args[1] - print """ ********************************************************************** ** ** @@ -213,7 +224,7 @@ def checkFile(name, comment): #locale.setlocale(locale.LC_ALL, 'fr_FR') time_start = datetime.datetime.now() -#print now,"(",now.strftime("%A %B %d %Y, %H:%M"),")" + print "Start time:", time_start @@ -248,7 +259,7 @@ def check_ffversion(reduced): allff_specs = { 'SCORPION': {'ff_file': 'scorpion.par', 'ff_class': ScorpionForceField, - 'minimizer_class': ScorpionLbfgs + 'minimizer_class': Lbfgs }, 'ATTRACT1': {'ff_file': 'aminon.par', @@ -287,22 +298,30 @@ def check_ffversion(reduced): - #load receptor and ligand: -rec=Rigidbody(options.receptor_name) -lig=Rigidbody(options.ligand_name) -rec=AttractRigidbody(rec) -lig=AttractRigidbody(lig) +if options.regions: + rec=Mcoprigid(options.receptor_name) + core=AttractRigidbody(options.ligand_name) + lig=Mcoprigid() + lig.setCore(core) +else: + rec=Rigidbody(options.receptor_name) + rec=AttractRigidbody(rec) + lig=Rigidbody(options.ligand_name) + lig=AttractRigidbody(lig) + print "Reading receptor (fixed): %s with %d particules" %( options.receptor_name, len(rec) ) print "Reading ligand (mobile): %s with %d particules" %( options.ligand_name, len(lig) ) -if (options.single and options.transnb): - parser.error("options -s and -t are mutually exclusive") +if (options.single and (options.transnb or options.rotnb)): + parser.error("options -s and -trans/rot are mutually exclusive") # save all minimization variables in trajectory file trjname = "minimization.trj" -if (options.single): +if (options.single or (options.transnb and options.rotnb)): ftraj = open(trjname, "w") + if(options.regions): + weight_variation = open("weight_variation.dat", "w") if (options.reffile): checkFile(options.reffile, "") @@ -315,6 +334,15 @@ def check_ffversion(reduced): else: Rmsd_alias = rmsdca +firstr = 0 +ff = open("translation.dat",'r') +for line in ff: + if "ATOM" in line: + spl = line.split() + firstr = int(spl[1])-1 + break +ff.close() + if (not options.single): #systematic docking with default translations and rotations # check for rotation.dat and translation.dat @@ -331,43 +359,44 @@ def check_ffversion(reduced): printFiles=True -# option -t used: define the selected translation +# option -trans used: define the selected translation transnb=0 if (options.transnb!=None): # check for rotation.dat and translation.dat checkFile("rotation.dat", "rotation file is required.") checkFile("translation.dat", "translation file is required.\nFormer users may rename translat.dat into translation.dat.") trans=Rigidbody("translation.dat") - - transnb=options.transnb - - if options.start1 is True: - transnb -= 1 - - co=trans.getCoords(transnb) - translations=[[transnb+1,co]] - - if transnb!= len(trans)-1: + co=trans.getCoords(options.transnb-1) + translations=[[options.transnb,co]] + transnb=options.transnb-1 + if transnb!=len(trans)-1: printFiles=False #don't append ligand, receptor, etc. unless this is the last translation point of the simulation - +if (options.rotnb!=None): + rotations = [rotations[options.rotnb-1]] # core attract algorithm for trans in translations: transnb+=1 - print "@@@@@@@ Translation nb %i @@@@@@@" %(transnb) + print "@@@@@@@ Translation nb %i @@@@@@@" %(transnb+firstr) rotnb=0 for rot in rotations: rotnb+=1 print "----- Rotation nb %i -----"%rotnb minimcounter=0 - ligand=AttractRigidbody(lig) + if options.regions: + ligand=Mcoprigid(lig) + receptor=Mcoprigid(rec) + else: + ligand=AttractRigidbody(lig) + receptor=AttractRigidbody(rec) center=ligand.FindCenter() ligand.Translate(Coord3D()-center) #set ligand center of mass to 0,0,0 ligand.AttractEulerRotate(surreal(rot[0]),surreal(rot[1]),surreal(rot[2])) ligand.Translate(trans[1]) + for minim in minimlist: minimcounter+=1 cutoff=math.sqrt(minim['squarecutoff']) @@ -377,41 +406,77 @@ def check_ffversion(reduced): #performs single minimization on receptor and ligand, given maxiter=niter and restraint constant rstk forcefield=ff_specs['ff_class'](ff_specs['ff_file'], surreal(cutoff) ) - rec.setTranslation(False) - rec.setRotation(False) + if options.regions: + mcopff = McopForceField(forcefield, surreal(cutoff)) + receptor.setTranslation(False) + receptor.setRotation(False) - forcefield.AddLigand(rec) - forcefield.AddLigand(ligand) + if options.regions: + mcopff.setLigand(ligand) + mcopff.setReceptor(receptor) + else: + forcefield.AddLigand(receptor) + forcefield.AddLigand(ligand) rstk=minim['rstk'] #restraint force #if rstk>0.0: #forcefield.SetRestraint(rstk) - lbfgs_minimizer=ff_specs['minimizer_class'](forcefield) - lbfgs_minimizer.minimize(niter) + + if options.regions: + lbfgs_minimizer=Lbfgs(mcopff) + lbfgs_minimizer.minimize(niter) + else : + lbfgs_minimizer=ff_specs['minimizer_class'](forcefield) + lbfgs_minimizer.minimize(niter) + X=lbfgs_minimizer.GetMinimizedVars() #optimized freedom variables after minimization - #TODO: test and use CenterToOrigin() ! - output=AttractRigidbody(ligand) + if options.regions: + output=Mcoprigid(ligand) + else: + output=AttractRigidbody(ligand) + center=output.FindCenter() output.Translate(Coord3D()-center) output.AttractEulerRotate(surreal(X[0]), surreal(X[1]), surreal(X[2])) output.Translate(Coord3D(surreal(X[3]),surreal(X[4]),surreal(X[5]))) output.Translate(center) - ligand=AttractRigidbody(output) - if (options.single): + if options.regions: + ligand=Mcoprigid(output) + else: + ligand=AttractRigidbody(output) + if (options.single or (options.transnb and options.rotnb)): ntraj=lbfgs_minimizer.GetNumberIter() for iteration in range(ntraj): traj = lbfgs_minimizer.GetMinimizedVarsAtIter(iteration) - for t in traj: + for t in traj[0:6]: ftraj.write("%f "%t) ftraj.write("\n") ftraj.write("~~~~~~~~~~~~~~\n") + if (options.regions): + savedWeights = mcopff.getSavedWeights() + for i, region in enumerate(savedWeights[0]): + weight_variation.write("Region %i\t"%(i+1)) + for k, copy in enumerate(region): + weight_variation.write("Copy %i\t"%(k+1)) + weight_variation.write("\n\n") + for iteration in range(ntraj): + for region in savedWeights[iteration]: + weight_variation.write("\t\t") + for copy in region: + weight_variation.write("%f\t"%copy) + weight_variation.write("\n") + weight_variation.write("~~~~~~~~~~~~~~\n") + #computes RMSD if reference structure available if (options.reffile): - rms=Rmsd_alias(ref, output) + if(options.regions): + rms=Rmsd_alias(ref, output.getCore()) + else: + rms=Rmsd_alias(ref, output) else: rms="XXXX" @@ -420,9 +485,14 @@ def check_ffversion(reduced): #with the new ligand position forcefield=ff_specs['ff_class'](ff_specs['ff_file'], surreal(500)) print "%4s %6s %6s %13s %13s" %(" ","Trans", "Rot", "Ener", "RmsdCA_ref") - pl = AttractPairList(rec, ligand,surreal(500)) - print "%-4s %6d %6d %13.7f %13s" %("==", transnb, rotnb, forcefield.nonbon8(rec,ligand,pl), str(rms)) - output.PrintMatrix() + if options.regions: + print "%-4s %6d %6d %13.7f %13s" %("==", transnb + firstr, rotnb, mcopff.CalcEnergy(receptor,ligand,forcefield,500), str(rms)) + output.getCore().PrintMatrix() #getCore because PrintMatrix works on AttractRigidy and not Mcoprigid + printWeights(receptor.getWeights(), mcopff.getMcopE()) + else: + pl = AttractPairList(receptor, ligand,surreal(500)) + print "%-4s %6d %6d %13.7f %13s" %("==", transnb + firstr, rotnb, forcefield.nonbon8(receptor,ligand,pl), str(rms)) + output.PrintMatrix() #output compressed ligand and receptor: @@ -435,8 +505,10 @@ def check_ffversion(reduced): print compress_file("attract.inp") # close trajectory file for single minimization -if (options.single): +if (options.single or (options.transnb and options.rotnb)): ftraj.close() + if (options.regions): + weight_variation.close() print "Saved all minimization variables (translations/rotations) in %s" %(trjname) # print end and elapsed time diff --git a/PyAttract/reduce.py b/PyAttract/reduce.py index 7a0df5d..232a9f5 100755 --- a/PyAttract/reduce.py +++ b/PyAttract/reduce.py @@ -6,6 +6,14 @@ import os from optparse import OptionParser +# ---- Functions ---- # +def is_digit_array(s): + for digit in s: + if not digit.isdigit(): + return False + return True +# -- End of Functions -- # + usage = "%prog --ff force_field [--dna] [--cgopt] [--dgrid 1.5] [--allow_missing] atomic_file.pdb > reduced_file.red" version = "%prog 1.0" parser = OptionParser(usage) @@ -25,6 +33,9 @@ parser.add_option("--allow_missing", action="store_true", dest="warning",default=False, help="don't stop program if atoms are missing, only display a warning on stderr") +# --mcop option: create multicopy reduced protein +parser.add_option("--mcop", dest="regions", default=False, help="positions of the multicopy region(s) separated by ':' (eg. --mcop 12-21:45-48)\n If entire protein, enter ':' (eg. --mcop :") + (options, args) = parser.parse_args() @@ -38,6 +49,8 @@ print "thispath is: ", thispath # define other parameter files + +# Check force field option (--ff) cmd_options=[] if options.ffversion not in ["attract1","attract2","scorpion"]: sys.stderr.write ("Error: please choose one of the following CG force field: attract1, attract2 or scorpion\n") @@ -47,12 +60,46 @@ cmd_options.append("--allow_missing") else: sys.stderr.write("Warning: allow_missing option not supported by SCORPION force field\n") + + +# Check multicopy option (--mcop) +if options.regions: + if options.regions.endswith('.pdb'): + cmd_options.append('--mcop .pdb') + args.insert(0,options.regions) + elif options.regions == ':': + cmd_options.append("--mcop :") + else: + pos = options.regions.replace(':','-').split('-') + reg = options.regions.split(':') + num_positions = len(pos) + num_regions = len(reg) + + if num_positions != 2*num_regions or not is_digit_array(pos): + sys.stderr.write ("Error: the --mcop input format is incorrect. \n Separate start and end of a region with '-'. \n Separate different regions with ':' \n (Example for 2 regions : --mcop 12-21:45-48)\n Or enter the files containing core, and region copies.\n (Example for 2 regions : --mcop core.pdb region1.pdb region2.pdb)\n") + sys.exit(1) + + starts = [] + ends = [] + for i in xrange(0,len(reg)): + starts.append(min(map(int, reg[i].split('-')))) + ends.append(max(map(int, reg[i].split('-')))) + starts, ends = zip(*sorted(zip(starts,ends))) + + for i in xrange(0,len(reg)-1): + if ends[i] >= starts[i+1]: + sys.stderr.write("Error: the --mcop input has overlapping regions.\n") + sys.exit(1) + + cmd_options.append("--mcop %s"%options.regions) + +#Launch the selected force field version. if options.ffversion=="scorpion": if options.optimizecharges: cmd_options.append("--cgopt") if options.delgrid: cmd_options.append("--dgrid %f"%options.delgrid) - cmd_options.append(args[0]) + cmd_options.extend(args) prgname = os.path.join(thispath, "reduce_scorpion.py") os.system(prgname + " " + " ".join(cmd_options)) @@ -62,13 +109,13 @@ cmd_options.append("--dna") else: cmd_options.append("--prot") - cmd_options.append(args[0]) + cmd_options.extend(args) prgname = os.path.join(thispath, "reduce_attract1.py") os.system(prgname + " %s"%" ".join(cmd_options)) if options.ffversion=="attract2": - cmd_options.append(args[0]) + cmd_options.extend(args) prgname = os.path.join(thispath, "reduce_attract2.py") os.system(prgname + " %s"%" ".join(cmd_options)) diff --git a/PyAttract/reduce_attract1.py b/PyAttract/reduce_attract1.py index 6e9a04a..dbb36ee 100755 --- a/PyAttract/reduce_attract1.py +++ b/PyAttract/reduce_attract1.py @@ -31,6 +31,9 @@ parser.add_option("--allow_missing", action="store_true", dest="warning",default=False, help="don't stop program if atoms are missing, only display a warning on stderr") +# --mcop option: create multicopy reduced protein +parser.add_option("--mcop", dest="regions", default=False, help="positions of the multicopy region(s) separated by ';' (eg. 12-21;45-48)") + (options, args) = parser.parse_args() @@ -297,91 +300,133 @@ def Show(self): #========================================================== # load atomic pdb file into Rigidbody object #========================================================== -allAtom=Rigidbody(atomicName) -sys.stderr.write("Load atomic file %s with %d atoms \n" %(atomicName, len(allAtom))) - -#extract all 'atoms' objects -atomList=[] -for i in xrange(len(allAtom)): - atom = allAtom.CopyAtom(i) - # look for residue or base type conversion - resName = atom.residType - if resName in resConv.keys(): - atom.residType = resConv[resName] - # look for atom type conversion - atomTag = atom.residType + '-' + atom.atomType - if atomTag in atomConv.keys(): - atomName = atomConv[atomTag].split('-')[1] - atom.atomType = atomName - atomList.append(atom) - -#count residues -residueTagList=[] -coarseResList=[] -for atom in atomList: - resName = atom.residType - # create a unique identifier for every residue - # resTag is for instance "LEU-296-A" - resTag = resName + '-'+ str(atom.residId) + '-' + atom.chainId - if resTag not in residueTagList: - if resBeadAtomModel.has_key(resName): - residueTagList.append(resTag) - # add a pattern residue to the list of coarse residues for the protein - # beware of the hugly list copy: use copy.deepcopy() ! - coarseResList.append(copy.deepcopy(resBeadAtomModel[resName])) - else: - sys.stderr.write("WARNING: residue %s is unknown the residues <-> beads <-> atoms list !!\n" %(resName)) - sys.stderr.write(" : residue %s will not be reduced into coarse grain\n" %(resName)) -sys.stderr.write("Number of residues: %i\n" %(len(residueTagList))) -#========================================================== -# iterate through all atoms and residues to fill beads -#========================================================== -sys.stderr.write("Reading all atoms and filling beads:\n") -for atom in atomList: - #resTag is like "LEU-296-A" - resTag = atom.residType + '-' + str(atom.residId) + '-' + atom.chainId - if resTag in residueTagList: - id = residueTagList.index(resTag) - coarseResList[id].FillAtom(atom.atomType, atom.coords.x, atom.coords.y, atom.coords.z) -#========================================================== -# reduce beads -#========================================================== -coarsegrainPdb = "" # complete coarse grain (reduced) pdb file -atomCnt = 0 # atom counter -sys.stderr.write("Coarse graining:\n") -for i in range(len(residueTagList)): - tag = residueTagList[i].split('-') - resName = tag[0] - resId = int(tag[1]) - coarseRes = coarseResList[i].Reduce(resName, resId) - for bead in coarseRes: - coord = bead[0] - atomName = bead[1] - atomTypeId = bead[2] - if atomTypeId in beadChargeDic: - atomCharge = beadChargeDic[atomTypeId] - else: - sys.stderr.write("WARNING: cannot find charge of bead %s %2d in %s \n" %(atomName, atomTypeId, ffName)) - sys.stderr.write(" : set default charge to 0.0\n") - atomCharge = 0.0 - prop = Atomproperty() - prop.atomType = atomName - atomCnt += 1 - prop.atomId = atomCnt - prop.residId = resId - prop.residType = resName - prop.chainId = ' ' - extra = ('%5i%8.3f%2i%2i') %(atomTypeId,atomCharge,0,0) - prop.extra = extra - newAtom = Atom(prop, coord) - coarsegrainPdb += newAtom.ToPdbString() + "\n" - -#========================================================== -# output coarse grain (reduced) pdb file -#========================================================== -sys.stdout.write("HEADER ATTRACT1 REDUCED PDB FILE\n") -sys.stdout.write(coarsegrainPdb) -sys.stderr.write("Coarse grain (reduced) output") -sys.stderr.write(": %d beads \n" %(atomCnt)) +models = [] +if options.regions: +# Preparing pdb multicopy models + if options.regions.endswith('.pdb'): + core = Rigidbody(args[0]) + region_copies = [] + for i in xrange(1, len(args)): + region_copies.append(Mcop(args[i])) + models.append(core) + models.append(region_copies) + write_mcop('_mcop.pdb', core, region_copies) + else: + core, region_copies = multicopy(atomicName, options.regions) + write_mcop('_mcop.pdb', core, region_copies) + models.append(core) + models.append(region_copies) +else: + models.append(Rigidbody(atomicName)) + sys.stderr.write("Load atomic file %s with %d atoms \n" %(atomicName, len(models[0]))) + +def reduce_model(allAtom): + #extract all 'atoms' objects + atomList=[] + for i in xrange(len(allAtom)): + atom = allAtom.CopyAtom(i) + # look for residue or base type conversion + resName = atom.residType + if resName in resConv.keys(): + atom.residType = resConv[resName] + # look for atom type conversion + atomTag = atom.residType + '-' + atom.atomType + if atomTag in atomConv.keys(): + atomName = atomConv[atomTag].split('-')[1] + atom.atomType = atomName + atomList.append(atom) + + #count residues + residueTagList=[] + coarseResList=[] + for atom in atomList: + resName = atom.residType + # create a unique identifier for every residue + # resTag is for instance "LEU-296-A" + resTag = resName + '-'+ str(atom.residId) + '-' + atom.chainId + if resTag not in residueTagList: + if resBeadAtomModel.has_key(resName): + residueTagList.append(resTag) + # add a pattern residue to the list of coarse residues for the protein + # beware of the hugly list copy: use copy.deepcopy() ! + coarseResList.append(copy.deepcopy(resBeadAtomModel[resName])) + else: + sys.stderr.write("WARNING: residue %s is unknown the residues <-> beads <-> atoms list !!\n" %(resName)) + sys.stderr.write(" : residue %s will not be reduced into coarse grain\n" %(resName)) + sys.stderr.write("Number of residues: %i\n" %(len(residueTagList))) + + #========================================================== + # iterate through all atoms and residues to fill beads + #========================================================== + sys.stderr.write("Reading all atoms and filling beads:\n") + for atom in atomList: + #resTag is like "LEU-296-A" + resTag = atom.residType + '-' + str(atom.residId) + '-' + atom.chainId + if resTag in residueTagList: + id = residueTagList.index(resTag) + coarseResList[id].FillAtom(atom.atomType, atom.coords.x, atom.coords.y, atom.coords.z) + #========================================================== + # reduce beads + #========================================================== + coarsegrainPdb = "" # complete coarse grain (reduced) pdb file + atomCnt = 0 # atom counter + sys.stderr.write("Coarse graining:\n") + for i in range(len(residueTagList)): + tag = residueTagList[i].split('-') + resName = tag[0] + resId = int(tag[1]) + coarseRes = coarseResList[i].Reduce(resName, resId) + for bead in coarseRes: + coord = bead[0] + atomName = bead[1] + atomTypeId = bead[2] + if atomTypeId in beadChargeDic: + atomCharge = beadChargeDic[atomTypeId] + else: + sys.stderr.write("WARNING: cannot find charge of bead %s %2d in %s \n" %(atomName, atomTypeId, ffName)) + sys.stderr.write(" : set default charge to 0.0\n") + atomCharge = 0.0 + prop = Atomproperty() + prop.atomType = atomName + atomCnt += 1 + prop.atomId = atomCnt + prop.residId = resId + prop.residType = resName + prop.chainId = ' ' + extra = ('%5i%8.3f%2i%2i') %(atomTypeId,atomCharge,0,0) + prop.extra = extra + newAtom = Atom(prop, coord) + coarsegrainPdb += newAtom.ToPdbString() + "\n" + + #========================================================== + # output coarse grain (reduced) pdb file + #========================================================== + sys.stdout.write(coarsegrainPdb) + sys.stderr.write("Coarse grain (reduced) output") + sys.stderr.write(": %d beads \n" %(atomCnt)) + +# iterate through all models (important for multicopy option) and write reduced pdb +print "HEADER ATTRACT1 REDUCED PDB FILE" +# iterate through 'models' varialbe which contains max 2 elements. +# The fist element of the models list is a Rigidbody (the core in the mcop case, and the whole protein in the normal case) +# In the mcop case, there is a second element, which is the region_copies list (a list of Mcop objects). +for i in xrange(0,len(models)): + if i == 0 and options.regions != ':': + sys.stderr.write("\nCore:\n") + print 'MODEL 0' + reduce_model(models[i]) + print 'ENDMDL 0' + if i >= 1: + # iterate through number of regions + for j in xrange(0, len(models[1])): + # iterate through number of regions copies + for k in xrange(0,len(models[1][j])): + sys.stderr.write("\nRegion %i Copy %i:\n" %(j+1, k+1)) + print 'MODEL ' + str(j+1) + ' ' + str(k+1) + reduce_model(models[1][j][k]) + print 'ENDMDL ' + str(j+1) + ' ' + str(k+1) + sys.stderr.write("\n\nSummary:\n") + for j in xrange(0, len(models[1])): + sys.stderr.write("Region %i: %i copies\n" %(j+1,len(models[1][j]))) diff --git a/PyAttract/reduce_attract2.py b/PyAttract/reduce_attract2.py index 5124660..7e8adc7 100755 --- a/PyAttract/reduce_attract2.py +++ b/PyAttract/reduce_attract2.py @@ -17,6 +17,9 @@ parser.add_option("--allow_missing", dest="ignoremissing", action="store_true", default=False, help="ignore missing heavy atoms (which will result in missing beads)" ) +# --mcop option: create multicopy reduced protein +parser.add_option("--mcop", dest="regions", default=False, help="positions of the multicopy region(s) separated by ';' (eg. 12-21;45-48)") + (options, args) = parser.parse_args() @@ -138,77 +141,118 @@ def create(self): +models = [] +if options.regions: +# Preparing pdb multicopy models + if options.regions.endswith('.pdb'): + core = Rigidbody(args[0]) + region_copies = [] + for i in xrange(1, len(args)): + region_copies.append(Mcop(args[i])) + models.append(core) + models.append(region_copies) + write_mcop('new_mcop.pdb', core, region_copies) + else: + core, region_copies = multicopy(args[0], options.regions) + write_mcop('new_mcop.pdb', core, region_copies) + models.append(core) + models.append(region_copies) +else: + models.append(Rigidbody(sys.argv[1])) + sys.stderr.write("Number of atoms:%i \n" %len(models[0]) ) + +def reduce_model(allAtom): + #extract all 'atoms' objects + atoms=[] + for i in xrange(len(allAtom)): + atoms.append(allAtom.CopyAtom(i)) + + + + #count residues: + residuMap={} + residulist=[] + + # chain id for the reduced file: + outChainId = atoms[0].chainId + + for at in atoms: + #fix for incorrect pdb: append a chainId when it's missing + if at.chainId=='': + at.chainId = 'A' + residueIdentifier = at.residType + str(at.chainId) + str(at.residId) + #residueIdentifier is like "LEUA296" + residuMap.setdefault(residueIdentifier, []).append(at) + if residueIdentifier not in residulist: + residulist.append(residueIdentifier) + + sys.stderr.write("Number of residues: %i\n" %(len(residuMap))) + sys.stderr.write("Start atom of each residue:\n") + orderedresid=[residuMap[i] for i in residulist ] + startatoms=[lat[0].atomId for lat in orderedresid ] + out = "" + for statom in startatoms: + out+=(str(statom))+" " + sys.stderr.write(out+"\n") + + + #iterates through all the residues and create reduced beads: + + totAtoms=0 + + index = 0 + for residKey, atomList in zip(residulist,orderedresid): + residType=residKey[:3] + if (residType)=="HIE": residType="HIS" #fix for an amber output file + #print "key:", residKey + residNumber=int(residKey[4:]) + #print "reducing residue %s of type %s" % (residKey, residType) + correspList=beadCorresp[residType] + #print correspList + for correspUnit in correspList: + atomTypeName=correspUnit[0] + lstToReduce=correspUnit[1] + atomTypeNumber=correspUnit[2] + atomCharge=correspUnit[3] + beadcreator=BeadCreator(atomTypeName,atomTypeNumber, atomCharge, lstToReduce, outChainId ) + for atom in atomList: + beadcreator.submit(atom) + try: + bead = beadcreator.create() + except IncompleteBead: + sys.stderr.write("The bead %i of residue %s is incomplete. Please check your pdb!\n"\ + %(totAtoms+1,residKey) ) + if not options.ignoremissing: raise + totAtoms+=1 + #now we must modify the bead: change the residue type and set the "extra" field correctly + bead.residType = residType + extra = ('%5i'+'%8.3f'+'%2i'*2) %(atomTypeNumber,atomCharge,0, 0) + bead.extra = extra + bead.atomId = totAtoms + bead.residId = residNumber + print bead.ToPdbString() + +# iterate through all models (important for multicopy option) and write reduced pdb +print "HEADER ATTRACT2 REDUCED PDB FILE" +# iterate through 'models' varialbe which contains max 2 elements. +# The fist element of the models list is a Rigidbody (the core in the mcop case, and the whole protein in the normal case) +# In the mcop case, there is a second element, which is the region_copies list (a list of Mcop objects). +for i in xrange(0,len(models)): + if i == 0 and options.regions != ':': + sys.stderr.write("\nCore:\n") + print 'MODEL 0' + reduce_model(models[i]) + print 'ENDMDL 0' + if i >= 1: + # iterate through number of regions + for j in xrange(0, len(models[1])): + # iterate through number of regions copies + for k in xrange(0,len(models[1][j])): + sys.stderr.write("\nRegion %i Copy %i:\n" %(j+1, k+1)) + print 'MODEL ' + str(j+1) + ' ' + str(k+1) + reduce_model(models[1][j][k]) + print 'ENDMDL ' + str(j+1) + ' ' + str(k+1) + sys.stderr.write("\n\nSummary:\n") + for j in xrange(0, len(models[1])): + sys.stderr.write("Region %i: %i copies\n" %(j+1,len(models[1][j]))) -allAtom=Rigidbody(sys.argv[1]) -sys.stderr.write("Number of atoms:%i \n" %len(allAtom) ) - -#extract all 'atoms' objects -atoms=[] -for i in xrange(len(allAtom)): - atoms.append(allAtom.CopyAtom(i)) - - - -#count residues: -residuMap={} -residulist=[] - -# chain id for the reduced file: -outChainId = atoms[0].chainId - -for at in atoms: - #fix for incorrect pdb: append a chainId when it's missing - if at.chainId=='': - at.chainId = 'A' - residueIdentifier = at.residType + str(at.chainId) + str(at.residId) - #residueIdentifier is like "LEUA296" - residuMap.setdefault(residueIdentifier, []).append(at) - if residueIdentifier not in residulist: - residulist.append(residueIdentifier) - -sys.stderr.write("Number of residues: %i\n" %(len(residuMap))) -sys.stderr.write("Start atom of each residue:\n") -orderedresid=[residuMap[i] for i in residulist ] -startatoms=[lat[0].atomId for lat in orderedresid ] -out = "" -for statom in startatoms: - out+=(str(statom))+" " -sys.stderr.write(out+"\n") - - -#iterates through all the residues and create reduced beads: - -totAtoms=0 - -print "HEADER ATTRACT2 REDUCED PDB FILE" -index = 0 -for residKey, atomList in zip(residulist,orderedresid): - residType=residKey[:3] - if (residType)=="HIE": residType="HIS" #fix for an amber output file - #print "key:", residKey - residNumber=int(residKey[4:]) - #print "reducing residue %s of type %s" % (residKey, residType) - correspList=beadCorresp[residType] - #print correspList - for correspUnit in correspList: - atomTypeName=correspUnit[0] - lstToReduce=correspUnit[1] - atomTypeNumber=correspUnit[2] - atomCharge=correspUnit[3] - beadcreator=BeadCreator(atomTypeName,atomTypeNumber, atomCharge, lstToReduce, outChainId ) - for atom in atomList: - beadcreator.submit(atom) - try: - bead = beadcreator.create() - except IncompleteBead: - sys.stderr.write("The bead %i of residue %s is incomplete. Please check your pdb!\n"\ - %(totAtoms+1,residKey) ) - if not options.ignoremissing: raise - totAtoms+=1 - #now we must modify the bead: change the residue type and set the "extra" field correctly - bead.residType = residType - extra = ('%5i'+'%8.3f'+'%2i'*2) %(atomTypeNumber,atomCharge,0, 0) - bead.extra = extra - bead.atomId = totAtoms - bead.residId = residNumber - print bead.ToPdbString() diff --git a/PyAttract/reduce_scorpion.py b/PyAttract/reduce_scorpion.py index 9046f87..e689605 100755 --- a/PyAttract/reduce_scorpion.py +++ b/PyAttract/reduce_scorpion.py @@ -3,15 +3,20 @@ from optparse import OptionParser +from cgopt import optimize parser = OptionParser() + +# --cgopt option: choice of charge optimization with SCORPION parser.add_option("--cgopt", dest="optimizecharges", action="store_true", default=False, help="optimize SCORPION coarse grained charges") parser.add_option("--dgrid", type="float",dest="delgrid", default=1.5, help="grid spacing (A) for charge optimization (default is 1.5), works only with -cgopt option") -(options, args) = parser.parse_args() +# --mcop option: create multicopy reduced protein +parser.add_option("--mcop", dest="regions", default=False, help="positions of the multicopy region(s) separated by ';' (eg. 12-21;45-48)") +(options, args) = parser.parse_args() import sys import copy @@ -153,171 +158,212 @@ def create(self): descriptions.append(beadDescription) beadCorresp[residname] = descriptions +models = [] +if options.regions: +# Preparing pdb multicopy models + if options.regions.endswith('.pdb'): + core = Rigidbody(args[0]) + region_copies = [] + for i in xrange(1, len(args)): + region_copies.append(Mcop(args[i])) + models.append(core) + models.append(region_copies) + write_mcop('new_mcop.pdb', core, region_copies) + else: + core, region_copies = multicopy(args[0], options.regions) + write_mcop('_mcop.pdb', core, region_copies) + models.append(core) + models.append(region_copies) +else: + models.append(Rigidbody(args[0])) + + + + +def reduce_model(allAtom): + newallAtom = [] + for i in xrange(len(allAtom)): + atom = allAtom.CopyAtom(i) + if atom.chainId == '': atom.chainId = ' ' + if atom.atomType[0] != 'H' and atom.atomType != 'OXT' and atom.atomType!= 'OT2': + newallAtom.append(atom) + allAtom = Rigidbody() + for at in newallAtom: + allAtom.AddAtom(at) + + + + sys.stderr.write("Number of atoms: %d\n" %(len(allAtom) )) + + #extract all 'atoms' objects + atoms=[] + for i in xrange(len(allAtom) ): + atoms.append(allAtom.CopyAtom(i)) + + + + #count residues: + residuMap={} + residulist=[] + for at in atoms: + residueIdentifier = at.residType + str(at.chainId) + str(at.residId) + #residueIdentifier is like "LEUA296" + residuMap.setdefault(residueIdentifier, []).append(at) + if residueIdentifier not in residulist: + residulist.append(residueIdentifier) + + sys.stderr.write("Number of residues: %i\n" %(len(residuMap))) + sys.stderr.write("Start atom of each residue:\n") + orderedresid=[residuMap[i] for i in residulist ] + startatoms=[lat[0].atomId for lat in orderedresid ] + out = "" + for statom in startatoms: + out+=(str(statom))+" " + sys.stderr.write(out+"\n") + + + #iterates through all the residues and create reduced beads: + + totAtoms=0 + + protein = [] + + for residKey, atomList in zip(residulist,orderedresid): + residType=residKey[:3] + if (residType)=="HIE": residType="HIS" #fix for an amber output file + residNumber=int(residKey[4:]) + correspList=beadCorresp[residType] + + for correspUnit in correspList: + atomTypeName=correspUnit[0] + lstToReduce=correspUnit[1] + atomTypeNumber=correspUnit[2] + atomCharge=correspUnit[3] + beadcreator=BeadCreator(atomTypeName,atomTypeNumber, atomCharge, lstToReduce) + for atom in atomList: + beadcreator.submit(atom) + try: + bead = beadcreator.create() + except IncompleteBead: + sys.stderr.write("The bead %i of residue %s is incomplete. Please check your pdb!\n"\ + %(totAtoms+1,residKey) ) + sys.exit(1) + totAtoms+=1 + #now we must modify the bead: change the residue type and set the "extra" field correctly + bead.residType = residType + extra = ('%5i'+'%8.3f'+'%2i'*2) %(atomTypeNumber,atomCharge,0, 0) + bead.extra = extra + bead.atomId = totAtoms + bead.residId = residNumber + protein.append(bead) + + + charge = [] + radius = [] + cx = [] + cy = [] + cz = [] + + cgch = [] + cgr = [] + cgx = [] + cgy = [] + cgz = [] + + + + for i in range(len(allAtom)): + atom = allAtom.CopyAtom(i) + residu_type= atom.residType + atomtype = atom.atomType + if residu_type =="ILE" and atomtype == "CD": + atomtype = "CD1" + key = "%s:%s"%(residu_type, atomtype) + radius.append( allAtomRadius[key] ) + charge.append ( allAtomCharges[key] ) + cx.append( atom.coords.x) + cy.append( atom.coords.y) + cz.append( atom.coords.z) + + + + for i, atom in enumerate(protein): + + cgch.append( atom.atomCharge ) + + residu_type= atom.residType + atomtype = atom.atomType + if residu_type =="ILE" and atomtype == "CD": + atomtype = "CD1" + key = "%s:%s"%(residu_type, atomtype) - -allAtom=Rigidbody(args[0]) - -newallAtom = [] -for i in xrange(len(allAtom)): - atom = allAtom.CopyAtom(i) - if atom.chainId == '': atom.chainId = ' ' - if atom.atomType[0] != 'H' and atom.atomType != 'OXT' and atom.atomType!= 'OT2': - newallAtom.append(atom) -allAtom = Rigidbody() -for at in newallAtom: - allAtom.AddAtom(at) - - - -sys.stderr.write("Number of atoms: %d\n" %(len(allAtom) )) - -#extract all 'atoms' objects -atoms=[] -for i in xrange(len(allAtom) ): - atoms.append(allAtom.CopyAtom(i)) - - - -#count residues: -residuMap={} -residulist=[] -for at in atoms: - residueIdentifier = at.residType + str(at.chainId) + str(at.residId) - #residueIdentifier is like "LEUA296" - residuMap.setdefault(residueIdentifier, []).append(at) - if residueIdentifier not in residulist: - residulist.append(residueIdentifier) - -sys.stderr.write("Number of residues: %i\n" %(len(residuMap))) -sys.stderr.write("Start atom of each residue:\n") -orderedresid=[residuMap[i] for i in residulist ] -startatoms=[lat[0].atomId for lat in orderedresid ] -out = "" -for statom in startatoms: - out+=(str(statom))+" " -sys.stderr.write(out+"\n") + cgr.append( beadRadius[key] ) + cgx.append( atom.coords.x) + cgy.append( atom.coords.y) + cgz.append( atom.coords.z) -#iterates through all the residues and create reduced beads: -totAtoms=0 -protein = [] -for residKey, atomList in zip(residulist,orderedresid): - residType=residKey[:3] - if (residType)=="HIE": residType="HIS" #fix for an amber output file - residNumber=int(residKey[4:]) - correspList=beadCorresp[residType] - for correspUnit in correspList: - atomTypeName=correspUnit[0] - lstToReduce=correspUnit[1] - atomTypeNumber=correspUnit[2] - atomCharge=correspUnit[3] - beadcreator=BeadCreator(atomTypeName,atomTypeNumber, atomCharge, lstToReduce) - for atom in atomList: - beadcreator.submit(atom) - try: - bead = beadcreator.create() - except IncompleteBead: - sys.stderr.write("The bead %i of residue %s is incomplete. Please check your pdb!\n"\ - %(totAtoms+1,residKey) ) - sys.exit(1) - totAtoms+=1 - #now we must modify the bead: change the residue type and set the "extra" field correctly - bead.residType = residType - extra = ('%5i'+'%8.3f'+'%2i'*2) %(atomTypeNumber,atomCharge,0, 0) - bead.extra = extra - bead.atomId = totAtoms - bead.residId = residNumber - protein.append(bead) + first = False + last = False + for i, at in enumerate(protein): + if at.atomType == 'CA': + cgch[i] += 1 + break + protein.reverse() + cgch.reverse() + for i,at in enumerate(protein): + if at.atomType == 'CA': + cgch[i] -= 1 + break + + protein.reverse() + cgch.reverse() + if options.optimizecharges: + optimized = optimize(len(allAtom), charge, radius, cx, cy, cz, len(protein), cgch, cgr, cgx, cgy, cgz, options.delgrid ) + else: + optimized = cgch -from cgopt import optimize + for i, bead in enumerate(protein): + #ugly hack to set correct charges values due to a bug atom to pdb conversion + extra = bead.extra + atomTypeNumber = int(extra.split()[0]) + bead.extra = ('%5i'+'%8.3f'+'%2i'*2) %(atomTypeNumber,optimized[i],0, 0) + + print bead.ToPdbString() -charge = [] -radius = [] -cx = [] -cy = [] -cz = [] - -cgch = [] -cgr = [] -cgx = [] -cgy = [] -cgz = [] - - - -for i in range(len(allAtom)): - atom = allAtom.CopyAtom(i) - residu_type= atom.residType - atomtype = atom.atomType - if residu_type =="ILE" and atomtype == "CD": - atomtype = "CD1" - key = "%s:%s"%(residu_type, atomtype) - radius.append( allAtomRadius[key] ) - charge.append ( allAtomCharges[key] ) - cx.append( atom.coords.x) - cy.append( atom.coords.y) - cz.append( atom.coords.z) - - - -for i, atom in enumerate(protein): - - cgch.append( atom.atomCharge ) - - residu_type= atom.residType - atomtype = atom.atomType - if residu_type =="ILE" and atomtype == "CD": - atomtype = "CD1" - key = "%s:%s"%(residu_type, atomtype) - - cgr.append( beadRadius[key] ) - cgx.append( atom.coords.x) - cgy.append( atom.coords.y) - cgz.append( atom.coords.z) - - - - - - -first = False -last = False -for i, at in enumerate(protein): - if at.atomType == 'CA': - cgch[i] += 1 - break -protein.reverse() -cgch.reverse() - -for i,at in enumerate(protein): - if at.atomType == 'CA': - cgch[i] -= 1 - break - -protein.reverse() -cgch.reverse() -if options.optimizecharges: - optimized = optimize(len(allAtom), charge, radius, cx, cy, cz, len(protein), cgch, cgr, cgx, cgy, cgz, options.delgrid ) -else: - optimized = cgch - +# iterate through all models (important for multicopy option) and write reduced pdb print "HEADER SCORPION REDUCED PDB FILE" - -for i, bead in enumerate(protein): - #ugly hack to set correct charges values due to a bug atom to pdb conversion - extra = bead.extra - atomTypeNumber = int(extra.split()[0]) - bead.extra = ('%5i'+'%8.3f'+'%2i'*2) %(atomTypeNumber,optimized[i],0, 0) - - print bead.ToPdbString() - +# iterate through 'models' variable which contains max 2 elements. +# The fist element of the models list is a Rigidbody (the core in the mcop case, and the whole protein in the normal case) +# In the mcop case, there is a second element, which is the region_copies list (a list of Mcop objects). +for i in xrange(0,len(models)): + if i == 0 and options.regions != ':': + sys.stderr.write("\nCore:\n") + print 'MODEL 0' + reduce_model(models[i]) + print 'ENDMDL 0' + if i >= 1: + # iterate through number of regions + for j in xrange(0, len(models[1])): + # iterate through number of regions copies + for k in xrange(0,len(models[1][j])): + sys.stderr.write("\nRegion %i Copy %i:\n" %(j+1, k+1)) + print 'MODEL ' + str(j+1) + ' ' + str(k+1) + reduce_model(models[1][j][k]) + print 'ENDMDL ' + str(j+1) + ' ' + str(k+1) + sys.stderr.write("\n\nSummary:\n") + for j in xrange(0, len(models[1])): + sys.stderr.write("Region %i: %i copies\n" %(j+1,len(models[1][j]))) + + + diff --git a/bindings/_ptools.pyx b/bindings/_ptools.pyx index d9a86e0..2d56744 100644 --- a/bindings/_ptools.pyx +++ b/bindings/_ptools.pyx @@ -10,8 +10,10 @@ include "atomselection.pyx" include "attractforcefield.pyx" include "scorpionforcefield.pyx" include "lbfgs.pyx" -include "scorpion_lbfgs.pyx" +#include "scorpion_lbfgs.pyx" include "rmsd.pyx" include "version.pyx" include "coordsarray.pyx" include "superpose.pyx" +include "mcopff.pyx" +include "multicopy.py" diff --git a/bindings/atomselection.pyx b/bindings/atomselection.pyx index 63f82d6..4b654e6 100644 --- a/bindings/atomselection.pyx +++ b/bindings/atomselection.pyx @@ -6,7 +6,7 @@ cdef extern from "atomselection.h" namespace "PTools": cdef cppclass CppAtomSelection "PTools::AtomSelection": CppAtomSelection() CppAtomSelection(CppAtomSelection&) - CppAtomSelection(CppRigidbody) + CppAtomSelection(CppRigidbody&) unsigned int Size() void SetRigid(CppRigidbody&) @@ -38,19 +38,20 @@ cdef class AtomSelection: self.thisptr = new CppAtomSelection() return - if isinstance(AtomSelection, arg): + elif isinstance(arg, AtomSelection): atsel = arg atselptr = atsel.thisptr self.thisptr = new CppAtomSelection(deref(atselptr)) return - if isinstance(Rigidbody, arg): + elif isinstance(arg, Rigidbody): rig = arg rigptr = rig.thisptr self.pyRigid = arg self.thisptr = new CppAtomSelection(deref(rigptr)) - raise RuntimeError("cannot reach here") + else: + raise RuntimeError("cannot reach here") def __dealloc__(self): if self.thisptr: @@ -95,4 +96,4 @@ cdef class AtomSelection: def SetRigid(self, Rigidbody r): self.pyRigid = r # to increase the refcount of r, preventing bad things if r is destroyed - self.thisptr.SetRigid(deref(r.thisptr)) \ No newline at end of file + self.thisptr.SetRigid(deref(r.thisptr)) diff --git a/bindings/attractforcefield.pyx b/bindings/attractforcefield.pyx index a600ccc..dd28d1f 100644 --- a/bindings/attractforcefield.pyx +++ b/bindings/attractforcefield.pyx @@ -1,6 +1,7 @@ from cython.operator cimport dereference as deref from libcpp.string cimport string + cdef extern from "" namespace "std": cdef cppclass vector[T]: cppclass iterator: @@ -37,11 +38,13 @@ cdef extern from "attractforcefield.h" namespace "PTools": cdef cppclass CppAttractForceField2 "PTools::AttractForceField2" (CppBaseAttractForceField) : CppAttractForceField2(string&, double) - -cdef class BaseAttractForceField: - cdef CppBaseAttractForceField * thisptr - - + +cdef class ForceField: + #cdef int junk_variable; # to avoid having empty cdef class + cdef CppForceField* thisptr + +cdef class BaseAttractForceField(ForceField): + cdef int junk_variable; # to avoid having empty cdef class cdef class AttractForceField2(BaseAttractForceField): @@ -49,6 +52,10 @@ cdef class AttractForceField2(BaseAttractForceField): #cdef CppAttractForceField2* thisptr def __cinit__(self, filename, cutoff): + # deallocate + del self.thisptr + self.thisptr = 0 + cdef char* c_filename cdef string * cppname @@ -63,23 +70,28 @@ cdef class AttractForceField2(BaseAttractForceField): def AddLigand(self, AttractRigidbody rig): self.rigidlist.append(rig) - self.thisptr.AddLigand(deref(rig.thisptr)) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + cpp_ptr.AddLigand(deref(rig.thisptr)) def Function(self, vec): cdef vector[double] v for el in vec: v.push_back(el) - return self.thisptr.Function(v) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.Function(v) def getVdw(self): - return self.thisptr.getVdw() + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.getVdw() def getCoulomb(self): - return self.thisptr.getCoulomb() + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.getCoulomb() def nonbon8(self, AttractRigidbody rec, AttractRigidbody lig, AttractPairList pl, verbose=False): - return self.thisptr.nonbon8(deref(rec.thisptr), deref(lig.thisptr), deref(pl.thisptr), verbose) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.nonbon8(deref(rec.thisptr), deref(lig.thisptr), deref(pl.thisptr), verbose) cdef extern from "attractforcefield.h" namespace "PTools": @@ -98,33 +110,43 @@ cdef class AttractForceField1(BaseAttractForceField): def __cinit__(self, filename, cutoff): + # deallocate + del self.thisptr + self.thisptr = 0 + cdef char* c_filename cdef string * cppname c_filename = filename cppname = new string(c_filename) - self.thisptr = new CppAttractForceField1(deref(cppname), cutoff) + self.thisptr = new CppAttractForceField1(deref(cppname), cutoff) del cppname def __dealloc__(self): del self.thisptr def AddLigand(self, AttractRigidbody rig): - self.thisptr.AddLigand(deref(rig.thisptr)) + self.rigidlist.append(rig) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + cpp_ptr.AddLigand(deref(rig.thisptr)) def Function(self, vec): cdef vector[double] v for el in vec: v.push_back(el) - return self.thisptr.Function(v) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.Function(v) def getVdw(self): - return self.thisptr.getVdw() + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.getVdw() def getCoulomb(self): - return self.thisptr.getCoulomb() + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.getCoulomb() def nonbon8(self, AttractRigidbody rec, AttractRigidbody lig, AttractPairList pl, verbose=False): - return self.thisptr.nonbon8(deref(rec.thisptr), deref(lig.thisptr), deref(pl.thisptr), verbose) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.nonbon8(deref(rec.thisptr), deref(lig.thisptr), deref(pl.thisptr), verbose) diff --git a/bindings/attractrigidbody.pyx b/bindings/attractrigidbody.pyx index 13942e1..3ebfdb0 100644 --- a/bindings/attractrigidbody.pyx +++ b/bindings/attractrigidbody.pyx @@ -46,9 +46,8 @@ cdef CppAttractRigidbody* _getAttractRigidbody_from_py_name(pyname): cdef class AttractRigidbody (Rigidbody) : - # cdef CppAttractRigidbody* thisptr - def __cinit__(self, arg): + def __cinit__(self, arg=''): # first deallocate the previously allocated Rigidbody del self.thisptr @@ -60,8 +59,12 @@ cdef class AttractRigidbody (Rigidbody) : self.thisptr = new CppAttractRigidbody(deref(rigidbodyptr)) return elif isinstance(arg, str): - self.thisptr = _getAttractRigidbody_from_py_name(arg) - return + if arg == '': + self.thisptr = new CppAttractRigidbody() + return + else: + self.thisptr = _getAttractRigidbody_from_py_name(arg) + return elif isinstance(arg, AttractRigidbody): oldrigidbody = arg oldrigidbody_ptr = oldrigidbody.thisptr @@ -145,4 +148,4 @@ cdef class AttractRigidbody (Rigidbody) : del ret.thisptr cdef CppAtomSelection new_sel = self.thisptr.CA() ret.thisptr = new CppAtomSelection(new_sel) - return ret \ No newline at end of file + return ret diff --git a/bindings/lbfgs.pyx b/bindings/lbfgs.pyx index fee4def..c486f67 100644 --- a/bindings/lbfgs.pyx +++ b/bindings/lbfgs.pyx @@ -3,11 +3,14 @@ from libcpp.string cimport string cdef extern from "lbfgs_interface.h" namespace "PTools": cdef cppclass CppLbfgs "PTools::Lbfgs": - CppLbfgs(CppBaseAttractForceField&) + CppLbfgs(CppForceField&) void minimize(int) vector[double] GetMinimizedVars() vector[double] GetMinimizedVarsAtIter(int) int GetNumberIter() + #void denormalize_weights() + #void normalize_weights() + vector[vector[double]] getWeights() cdef class Lbfgs: @@ -15,8 +18,8 @@ cdef class Lbfgs: def __cinit__(self, forcefield): - ff = forcefield - cdef CppBaseAttractForceField * ffptr = ff.thisptr + ff = forcefield + cdef CppForceField * ffptr = ff.thisptr self.thisptr = new CppLbfgs(deref(ffptr)) def __dealloc__(self): @@ -40,4 +43,13 @@ cdef class Lbfgs: out = [] for i in xrange(vars.size()): out.append(vars[i]) - return out \ No newline at end of file + return out + + #def denormalize_weights(self): + # self.thisptr.denormalize_weights() + + #def normalize_weights(self): + # self.thisptr.normalize_weights() + + def getWeights(self): + self.thisptr.getWeights() \ No newline at end of file diff --git a/bindings/mcopff.pyx b/bindings/mcopff.pyx new file mode 100644 index 0000000..2e4977c --- /dev/null +++ b/bindings/mcopff.pyx @@ -0,0 +1,333 @@ +from cython.operator cimport dereference as deref +from libcpp.string cimport string +from cpython cimport bool + +cdef extern from "" namespace "std": + cdef cppclass stringstream: + stringstream() + string str() + void str(string &s) + +cdef extern from "mcopff.h" namespace "PTools": + cdef cppclass CppMcop "PTools::Mcop": + CppMcop() + CppMcop(string&) except+ + CppMcop(CppMcop&) + CppRigidbody getCopy(unsigned int) + CppRigidbody operator[](unsigned int) + unsigned int size() + void addCopy(CppRigidbody&) + cdef cppclass CppAttractMcop "PTools::AttractMcop": + CppAttractMcop() + CppAttractMcop(string&) except+ + CppAttractMcop(CppMcop&) + CppAttractMcop(CppAttractMcop&) + CppAttractRigidbody getCopy(unsigned int) + CppAttractRigidbody operator[](unsigned int) + unsigned int size() + void addCopy(CppAttractRigidbody&) + cdef cppclass CppMcoprigid "PTools::Mcoprigid": + CppMcoprigid() + CppMcoprigid(string&) except+ + CppMcoprigid(CppMcoprigid&) + CppAttractMcop getRegion(unsigned int) + CppAttractRigidbody getCore() + unsigned int size() + void setCore(CppAttractRigidbody&) + vector[vector[double]] getWeights() + Coord3D FindCenter() + void setTranslation(int) + void setRotation(int) + bool checkTranslation() + bool checkRotation() + void Translate(CppCoord3D&) + void AttractEulerRotate(double, double, double) + + cdef cppclass CppMcopForceField "PTools::McopForceField": + CppMcopForceField(CppBaseAttractForceField&, double) + vector[vector[double]] getWeights() + vector[vector[double]] getMcopE() + void setReceptor(CppMcoprigid&) + void setLigand(CppMcoprigid&) + double CalcEnergy(CppMcoprigid&, CppMcoprigid&, CppBaseAttractForceField&, double) + vector[vector[vector[double]]] getSavedWeights() + +cdef class Mcop: + cdef CppMcop* thisptr + + def __cinit__(self, filename=''): + + cdef CppMcop* oldptr + cdef Mcop oldmcop + cdef char* name + cdef string *cppname + cdef CppMcop* newmcop + + if isinstance(filename, str): + if filename == '': + self.thisptr = new CppMcop() + else: + # there is a filename, loading the pdb file + name = filename + cppname = new string(name) + newmcop = new CppMcop(deref(cppname)) + del cppname + self.thisptr = newmcop + + elif isinstance(filename, Mcop): + oldmcop = filename + oldptr = (oldmcop.thisptr) + self.thisptr = new CppMcop(deref(oldptr)) + if not self.thisptr: + print "FATAL: this should never happen" + + else: + raise RuntimeError("invalid argument in Mcop()") + + def __dealloc__(self): + if self.thisptr: + del self.thisptr + self.thisptr = 0 + + def getCopy(self, unsigned int i): + cdef CppRigidbody cpp_copy = self.thisptr.getCopy(i) + cdef Rigidbody cop = Rigidbody() + cdef CppRigidbody* new_cop = new CppRigidbody(cpp_copy) + del cop.thisptr + cop.thisptr = new_cop + return cop + + def __getitem__(self,unsigned int i): + cdef CppRigidbody cpp_copy = self.thisptr.getCopy(i) + cdef Rigidbody cop = Rigidbody() + cdef CppRigidbody* new_cop = new CppRigidbody(cpp_copy) + del cop.thisptr + cop.thisptr = new_cop + return cop + + def __len__(self): + return self.thisptr.size() + + def addCopy(self, Rigidbody rigid): + self.thisptr.addCopy(deref(rigid.thisptr)) + + +cdef class AttractMcop (Mcop): + cdef CppAttractMcop* thisptr2 + + def __cinit__(self, filename=''): + + # first deallocate the previously allocated Mcop + del self.thisptr + self.thisptr = new CppMcop() + + cdef CppAttractMcop* oldptr + cdef CppMcop* oldmcop_ptr + cdef AttractMcop oldattractmcop + cdef Mcop oldmcop + cdef char* name + cdef string *cppname + cdef CppAttractMcop* newattractmcop + + if isinstance(filename, str): + if filename == '': + self.thisptr2 = new CppAttractMcop() + else: + # there is a filename, laoding the pdb file + + name = filename + cppname = new string(name) + newattractmcop = new CppAttractMcop(deref(cppname)) + del cppname + self.thisptr2 = newattractmcop + + elif isinstance(filename, AttractMcop): + oldattractmcop = filename + oldptr = (oldattractmcop.thisptr2) + self.thisptr2 = new CppAttractMcop(deref(oldptr)) + if not self.thisptr2: + print "FATAL: this should never happen" + + elif isinstance(filename, Mcop): + + oldmcop = filename + oldmcop_ptr = (oldmcop.thisptr) + self.thisptr2 = new CppAttractMcop(deref(oldmcop_ptr)) + if not self.thisptr2: + print "FATAL; this should never happen" + + else: + raise RuntimeError("invalid argument in AttractMcop()") + + def __dealloc__(self): + if self.thisptr: + del self.thisptr + self.thisptr = 0 + if self.thisptr2: + del self.thisptr2 + self.thisptr2 = 0 + + def getCopy(self, unsigned int i): + cdef CppAttractRigidbody cpp_copy = self.thisptr2.getCopy(i) + cdef AttractRigidbody cop = AttractRigidbody() + cdef CppAttractRigidbody* new_cop = new CppAttractRigidbody(cpp_copy) + del cop.thisptr + cop.thisptr = new_cop + return cop + + def __getitem__(self,unsigned int i): + cdef CppAttractRigidbody cpp_copy = self.thisptr2.getCopy(i) + cdef AttractRigidbody cop = AttractRigidbody() + cdef CppAttractRigidbody* new_cop = new CppAttractRigidbody(cpp_copy) + del cop.thisptr + cop.thisptr = new_cop + return cop + + def addCopy(self, AttractRigidbody attrigid): + self.thisptr2.addCopy(deref(attrigid.thisptr)) + + def __len__(self): + return self.thisptr2.size() + + +cdef class Mcoprigid: + cdef CppMcoprigid* thisptr + + def __cinit__(self, filename='', arg2=''): + + cdef CppMcoprigid* oldptr + cdef Mcoprigid oldmcoprigid + cdef char* name + cdef string *cppname + cdef CppMcoprigid* newmcoprigid + + if isinstance(filename, str): + if filename == '': + self.thisptr = new CppMcoprigid() + else: + # there is a filename, loading the pdb file + name = filename + cppname = new string(name) + newmcoprigid = new CppMcoprigid(deref(cppname)) + del cppname + self.thisptr = newmcoprigid + + elif isinstance(filename, Mcoprigid): + oldmcoprigid = filename + oldptr = oldmcoprigid.thisptr + self.thisptr = new CppMcoprigid(deref(oldptr)) + if not self.thisptr: + print "FATAL: this should never happen" + + #if filename is AttractRigibody and arg2 is list of AttractMcop + #elif isinstance(arg2, AttractRigidbody) and isinstance(filename, (list, tuple)) and all(isinstance(elem, AttractMcop) for elem in arg2): + # TODO + + else: + raise RunetimeError("invalid argument in Mcoprigid()") + + def __dealloc__(self): + if self.thisptr: + del self.thisptr + self.thisptr = 0 + + def getCore(self): + cdef CppAttractRigidbody cpp_core = self.thisptr.getCore() + cdef AttractRigidbody core = AttractRigidbody() + cdef CppAttractRigidbody* new_core = new CppAttractRigidbody(cpp_core) + del core.thisptr + core.thisptr = new_core + return core + + def getRegion(self,i): + cdef CppAttractMcop cpp_region = self.thisptr.getRegion(i) + cdef AttractMcop region = AttractMcop() + cdef CppAttractMcop* new_region = new CppAttractMcop(cpp_region) + del region.thisptr + region.thisptr2 = new_region + return region + + def setCore(self, AttractRigidbody core): + self.thisptr.setCore(deref(core.thisptr)) + + def getWeights(self): + return self.thisptr.getWeights() + + + def FindCenter(self): + return self.getCore().FindCenter() + + def setTranslation(self, flag): + self.thisptr.setTranslation(flag) + + def setRotation(self, flag): + self.thisptr.setRotation(flag) + + def checkTranslation(self): + return self.thisptr.checkTranslation() + + def checkRotation(self): + return self.thisptr.checkRotation(); + + def AttractEulerRotate(self, double phi, double ssi, double rot): + self.thisptr.AttractEulerRotate(phi, ssi, rot) + + def Translate(self, Coord3D tr): + self.thisptr.Translate(deref(tr.thisptr)) + + def __len__(self): + return self.thisptr.size() + + +cdef class McopForceField(ForceField): + + def __cinit__(self, ff, cutoff): + # deallocate + del self.thisptr + self.thisptr = 0 + + cdef BaseAttractForceField old_ff + cdef CppBaseAttractForceField c_ff + cdef CppBaseAttractForceField* c_ff_ptr + + old_ff = ff + c_ff_ptr = old_ff.thisptr + self.thisptr = new CppMcopForceField(deref(c_ff_ptr), cutoff) + + + def __dealloc__(self): + if self.thisptr: + del self.thisptr + self.thisptr = 0 + + def getWeights(self): + cdef CppMcopForceField* cpp_ptr = self.thisptr + return cpp_ptr.getWeights() + + def getMcopE(self): + cdef CppMcopForceField* cpp_ptr = self.thisptr + return cpp_ptr.getMcopE() + + def setReceptor(self, Mcoprigid rec): + cdef CppMcopForceField* cpp_ptr = self.thisptr + cpp_ptr.setReceptor(deref(rec.thisptr)) + + def setLigand(self, Mcoprigid lig): + cdef CppMcopForceField* cpp_ptr = self.thisptr + cpp_ptr.setLigand(deref(lig.thisptr)) + + def CalcEnergy(self, Mcoprigid rec, Mcoprigid lig, BaseAttractForceField ff, double cutoff): + + cdef CppMcopForceField* cpp_ptr = self.thisptr + cdef CppBaseAttractForceField* cppff + + cppff = ff.thisptr + return cpp_ptr.CalcEnergy(deref(rec.thisptr), deref(lig.thisptr), deref(cppff), cutoff) + + def getSavedWeights(self): + + cdef CppMcopForceField* cpp_ptr = self.thisptr + return cpp_ptr.getSavedWeights() + + + \ No newline at end of file diff --git a/bindings/multicopy.py b/bindings/multicopy.py new file mode 100644 index 0000000..ef66f36 --- /dev/null +++ b/bindings/multicopy.py @@ -0,0 +1,265 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +from ptools import * +import sys + + +def multicopy(models_pdb, mcop_arg): + """ + Takes the filename of a pdb containing multiple models and the string of + the mcop argument indicating the region positions. + + Returns a single core of the protein (Rigidbody object) and all the copies + of the region(s) (list of Mcop objects where each Mcop object contains the + copies of a region which are Rigidbody) with respect to the mcop argument. + + ex: core, region_copies = ("NMRmodels.pdb", "23-44:66-79") + + The mcop argument uses the following format: + - If there are different regions, they are separated by ':' + - The start and end position of a region are separated by '-' + - To fetch copies of the whole protein, the mcop argument is ':' + + """ + starts, ends = convert_mcop_argument(mcop_arg) + models = Mcop(models_pdb) + cores, region_copies = extract_regions(models, starts, ends) + if mcop_arg == ':': + core = Rigidbody() + moved_copies = region_copies + else: + core = find_core(cores) + matrices = get_align_matrices(cores, core) + moved_copies = [] + for i in xrange(0,len(region_copies)): + moved_copies.append(apply_matrices(region_copies[i], matrices)) + return core, moved_copies + + +def convert_mcop_argument(arg): + """ + Takes the mcop arguement (string). + + Returns start and end positions (list of int) of the region(s). + + ex: starts, ends = convert_mcop_argument("23-44:66-79") + - starts variable contains [23, 66] + - ends variable contains [66, 79] + + Note : if the argument is ':', returns (-1,) and (-1,). This case is for + fetching copies of the whole protein. See extract_regions() + """ + # if whole protein + if arg == ':': + return ((-1,), (-1,)) + else: + reg = arg.split(':') + starts = [] + ends = [] + for i in xrange(0,len(reg)): + starts.append(min(map(int, reg[i].split('-')))) + ends.append(max(map(int, reg[i].split('-')))) + starts, ends = zip(*sorted(zip(starts,ends))) + return starts, ends + + +def extract_regions(models, starts, ends): + """ + Takes different models of a protein (Mcop object), start + position(s) of region(s) (list int), and end position(s) of + region(s) (list int). + + Returns the different cores of the protein (Mcop object) + and all the copies of the region(s) (list of Mcop objects where each Mcop + object contains the copies of a region) with respect to the start and + end positions. + + ex: cores, region_copies = extract_regions(models, starts, ends) + + Note: To fetch copies of the whole protein, starts and ends must both + contain -1 (int) at index 0. + + Note: a single core can then be selected from the returned multiple cores + thanks to the find_core() function. + """ + mcop = [] + cores = Mcop() + # loop through models + for i in xrange(0,len(models)): + reg = [] + model = models[i] + # selecting first region + reg.append(model.SelectResRange(starts[0], ends[0])) + fused_regions = reg[0] + if i == 0: + mcop.append(Mcop()) + #Creating and adding copy of region + copy = reg[0] + mcop[0].addCopy(copy.CreateRigid()) + for j in xrange(1, len(starts)): + #adding other region + if i == 0: + mcop.append(Mcop()) + #Creating and adding copy of region + reg.append(model.SelectResRange(starts[j], ends[j])) + copy = reg[j] + mcop[j].addCopy(copy.CreateRigid()) + fused_regions = copy | fused_regions + #Creating and Adding core (inverse of regions) + core = AtomSelection(model) + core = core & ~fused_regions + core = core.CreateRigid() + cores.addCopy(core) + # if whole protein + if starts[0] == -1 and ends[0] == -1: + temp = cores + cores = mcop[0] + mcop[0] = temp + return cores, mcop + + +def write_region_copies(filename, models, append=False): + """ + Writes a list of Mcop objects in a pdb file where each written model is + delimited by MODEL i j and ENDMDL i j. (i = region number, j = copy number) + + Takes an output filename (string), and a list of Mcop objects. + + Default arguments : + - append option (False if overide, True if append), default value = False + """ + if append: + f = open(filename, 'a') + else: + f = open(filename, 'w') + # iterate through number of regions + for i in xrange(0, len(models)): + # iterate through number of regions copies + for j in xrange(0,len(models[i])): + f.write('MODEL ' + str(i+1) + ' ' + str(j+1) + '\n') + f.write(str(models[i][j]) + '\n') + f.write('ENDMDL ' + str(i+1) + ' ' + str(j+1) + '\n') + f.close() + + +def write_single_model(filename, model, append=False, model_index=1): + """ + Writes a single Rigidbody object into a pdb file. The written model is + delimited by MODEL i and ENDMDL i (i = index value of the model). + + Takes an output filename (string), and a Rigidbody object. + + Default arguments : + - "append": False if overide file, True if append file, default value = False + - "model_index": index value of the written model (int), default value = 1 + """ + if append: + f = open(filename, 'a') + else: + f = open(filename, 'w') + f.write('MODEL ' + str(model_index) + '\n') + f.write(str(model) + '\n') + f.write('ENDMDL ' + str(model_index) + '\n') + f.close() + + +def write_mcop(filename, core, region_copies): + """ + Writes a multicopy pdb file. This type of file uses the following format : + The first model is the core of the protein and is delimited by MODEL 0 and + ENDMDL 0. + + The following models are the copies of the defined variable region(s). + These are delimited by MODEL i j and ENDMDL i j + (i = region number, j = copy number). + + Takes an output filename (string), the protein core (Rigidbody object), + and the region copies (list of Mcop objects). + """ + if len(core) > 0: + write_single_model(filename, core, model_index=0) + write_region_copies(filename, region_copies, append=True) + + +def rmsd_sums(cores): + """ + Does pairwise superpositions of proteins, and for each protein, calculates + the sum of the RMSD obtained from all of its superpositions. + + Takes a list of proteins, generally protein cores + (list of Rigidbody objects). + + Returns a list of the corresponding RMSD sums (list of int). + """ + rmsd_sums = [0 for x in xrange(0,len(cores))] + for i in xrange(0, len(cores)): + core1 = cores[i] + for j in xrange(i+1,len(cores)): + core2 = cores[j] + sup = superpose(core1.CA().CreateRigid(), core2.CA().CreateRigid()) + rmsd_sums[i] += sup.rmsd + rmsd_sums[j] += sup.rmsd + return rmsd_sums + + +def find_core(cores): + """ + From a list of proteins, returns the protein with the smallest RMSD sum + (with pairwise superpositions using the rmsd_sums function). This + is generally used to select the protein core from a list of potential cores. + + Takes of a list of Rigidbody objects. + + Returns a single Rigidbody object. + """ + sys.stderr.write("\nFinding core:\n") + sums = rmsd_sums(cores) + core_index = sums.index(min(sums)) + for i,sum in enumerate(sums): + sys.stderr.write("Core candidate %i: " %(i+1)) + sys.stderr.write("mean RMSD %f A \n" %(sum/(len(sums)-1))) + sys.stderr.write("\nSelected core candiate %i: mean RMSD %f A\n"\ + %(core_index+1, sums[core_index]/(len(sums)-1))) + return cores[core_index] + + +def get_align_matrices(cores, core): + """ + Finds the displacement matrices (translation and rotation) of + superpositions of a list of proteins (generally counterselected + protein cores) over one protein (generally the selected protein core). + These matrices are usefull to then move copies of other regions of the + protein (they generally correspond to, but are not present in, the + conterselected protein cores) to better position them onto the selected + protein core. + + Takes a Mcop objects, and a single Rigidbody objects. + Returns a list of Matrix objects corresponding to the entered list of + Rigidbody objects. + """ + matrices = [] + for i in xrange(0, len(cores)): + mobile_core = cores[i] + sup = superpose(core, mobile_core) + matrices.append(sup.matrix) + return matrices + + +def apply_matrices(region_copies, matrices): + """ + Modifies the coordinates of a list of Mcop objects according to a list of + Matrix objects. This function is usually used after obtaining the matrices + from the get_align_matrices function. + + Takes a series of region copies (list of Mcop objects) + Returns a series of moved region copies (list of Mcop ojects with + new cordinates) + + """ + moved_copies = Mcop() + for i in xrange(0, len(region_copies)): + cop = region_copies[i] + cop.ApplyMatrix(matrices[i]) + moved_copies.addCopy(cop) + return moved_copies + diff --git a/bindings/scorpion_lbfgs.pyx b/bindings/scorpion_lbfgs.pyx deleted file mode 100644 index 0ac2a22..0000000 --- a/bindings/scorpion_lbfgs.pyx +++ /dev/null @@ -1,43 +0,0 @@ -from cython.operator cimport dereference as dereference -from libcpp.string cimport string - -cdef extern from "lbfgs_interface.h" namespace "PTools": - cdef cppclass CppScorpionLbfgs "PTools::Lbfgs": - CppScorpionLbfgs(CppScorpionForceField&) - void minimize(int) - vector[double] GetMinimizedVars() - vector[double] GetMinimizedVarsAtIter(int) - int GetNumberIter() - -cdef class ScorpionLbfgs: - - cdef CppScorpionLbfgs * thisptr - - def __cinit__(self, forcefield): - - ff = forcefield - cdef CppScorpionForceField * ffptr = ff.thisptr - self.thisptr = new CppScorpionLbfgs(deref(ffptr)) - - def __dealloc__(self): - del self.thisptr - - def minimize(self, int maxiter): - self.thisptr.minimize(maxiter) - - def GetMinimizedVars(self): - cdef vector[double] vars = self.thisptr.GetMinimizedVars() - out = [] - for i in xrange(vars.size()): - out.append(vars[i]) - return out - - def GetNumberIter(self): - return self.thisptr.GetNumberIter() - - def GetMinimizedVarsAtIter(self, int iter): - cdef vector[double] vars = self.thisptr.GetMinimizedVarsAtIter(iter) - out = [] - for i in xrange(vars.size()): - out.append(vars[i]) - return out diff --git a/bindings/scorpionforcefield.pyx b/bindings/scorpionforcefield.pyx index 0be6502..d172234 100644 --- a/bindings/scorpionforcefield.pyx +++ b/bindings/scorpionforcefield.pyx @@ -2,6 +2,14 @@ from cython.operator cimport dereference as deref from libcpp.string cimport string +cdef extern from "forcefield.h" namespace "PTools": + cdef cppclass CppForceField "PTools::ForceField": + pass + +cdef extern from "attractforcefield.h" namespace "PTools": + cdef cppclass CppBaseAttractForceField "PTools::BaseAttractForceField": + pass + cdef extern from "scorpionforcefield.h" namespace "PTools": cdef cppclass CppScorpionForceField "PTools::ScorpionForceField": CppScorpionForceField(string&, double) @@ -12,38 +20,48 @@ cdef extern from "scorpionforcefield.h" namespace "PTools": double nonbon8(CppAttractRigidbody& , CppAttractRigidbody& , CppAttractPairList & , int) -cdef class ScorpionForceField: +cdef class ScorpionForceField (BaseAttractForceField): - cdef CppScorpionForceField* thisptr + #cdef CppScorpionForceField* thisptr def __cinit__(self, filename, cutoff): + + # deallocate + del self.thisptr + self.thisptr = 0 + cdef char* c_filename cdef string * cppname c_filename = filename cppname = new string(c_filename) - self.thisptr = new CppScorpionForceField(deref(cppname), cutoff) + self.thisptr = new CppScorpionForceField(deref(cppname), cutoff) del cppname def __dealloc__(self): del self.thisptr def AddLigand(self, AttractRigidbody rig): - self.thisptr.AddLigand(deref( rig.thisptr)) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + cpp_ptr.AddLigand(deref(rig.thisptr)) def Function(self, vec): cdef vector[double] v for el in vec: v.push_back(el) - return self.thisptr.Function(v) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.Function(v) def getVdw(self): - return self.thisptr.getVdw() + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.getVdw() def getCoulomb(self): - return self.thisptr.getCoulomb() + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.getCoulomb() def nonbon8(self, AttractRigidbody rec, AttractRigidbody lig, AttractPairList pl, verbose=False): - return self.thisptr.nonbon8(deref(rec.thisptr), deref(lig.thisptr), deref(pl.thisptr), verbose) + cdef CppBaseAttractForceField* cpp_ptr = self.thisptr + return cpp_ptr.nonbon8(deref(rec.thisptr), deref(lig.thisptr), deref(pl.thisptr), verbose) diff --git a/headers/attractforcefield.h b/headers/attractforcefield.h index 76c35a4..882b5aa 100644 --- a/headers/attractforcefield.h +++ b/headers/attractforcefield.h @@ -40,7 +40,6 @@ class BaseAttractForceField: public ForceField { std::vector forcesrec (rec.Size()); std::vector forceslig (lig.Size()); - dbl ener = nonbon8_forces(rec, lig, pairlist, forcesrec, forceslig, print); rec.addForces(forcesrec); lig.addForces(forceslig); diff --git a/headers/attractrigidbody.h b/headers/attractrigidbody.h index e014ae2..ed76c9b 100644 --- a/headers/attractrigidbody.h +++ b/headers/attractrigidbody.h @@ -57,7 +57,6 @@ class AttractRigidbody: public Rigidbody }; virtual bool isAtomActive(uint i) const { - uint atomtype = this->m_atomTypeNumber[i]; for(uint j=0; j& dummy); ///< set a list of ignored atom types @@ -87,6 +88,9 @@ class AttractRigidbody: public Rigidbody void updateActiveList(); + void AddAtom(const Atom& at); + void AddAtom(const Atomproperty& at, Coord3D co); + private: void init_(); diff --git a/headers/coordsarray.h b/headers/coordsarray.h index 22b0dae..fd28c03 100644 --- a/headers/coordsarray.h +++ b/headers/coordsarray.h @@ -21,7 +21,7 @@ #define _COORDS_ARRAY_H_ #include -#include +#include #include "coord3d.h" diff --git a/headers/forcefield.h b/headers/forcefield.h index f47408a..b37d1b2 100644 --- a/headers/forcefield.h +++ b/headers/forcefield.h @@ -42,6 +42,10 @@ class ForceField ///virtual destructor (Effective C++ - Scott Meyers - Item 7) virtual ~ForceField(){}; + + /// Saves weights into attribute. Used at each minimization iteration in order to follow + /// the weights' evolution during the docking. Used only for McopForcefield. + virtual void saveWeights(){}; } ; diff --git a/headers/lbfgs_interface.h b/headers/lbfgs_interface.h index a149275..d81ef76 100644 --- a/headers/lbfgs_interface.h +++ b/headers/lbfgs_interface.h @@ -27,6 +27,10 @@ class Lbfgs std::vector GetMinimizedVarsAtIter(uint iter); int GetNumberIter() {return m_opt->niter;} + //void denormalize_weights(); + //void normalize_weights(); + std::vector< std::vector > getWeights(); + diff --git a/headers/mcopff.h b/headers/mcopff.h index 4e5048a..6fa915b 100644 --- a/headers/mcopff.h +++ b/headers/mcopff.h @@ -9,48 +9,71 @@ #include "attractforcefield.h" #include "rigidbody.h" +#include "pdbio.h" namespace PTools{ // typedef std::vector ensemble ; - -class Region +class Mcop { private: -std::vector _copies; +std::vector _copies; public: -/* - void AttractEulerRotate(const dbl& phi, const dbl& ssi, const dbl& rot) - { - for (uint i=0; i<_copies.size(); i++) - { - _copies[i].AttractEulerRotate(phi, ssi, rot); - } - }; + //Contructor : no arguments + Mcop(){ } + //Contructor : takes filename + Mcop(std::string filename); + + + virtual void addCopy(const Rigidbody& cop){_copies.push_back(cop);}; - void Translate(const Coord3D& co) - { - for(uint i =0; i< _copies.size(); ++i) - { - _copies[i].Translate(co); - } - }*/ + + virtual size_t size() const {return _copies.size();}; + + virtual Rigidbody& operator[](uint i){return _copies[i];}; + virtual Rigidbody& getCopy(int i){return _copies[i];}; - void addCopy(const AttractRigidbody& cop){_copies.push_back(cop);}; + void ReadmcopPDB(const std::string name); + void ReadmcopPDB(std::istream& file, std::vector& protein); + void ReadModelsPDB(std::istream& file, std::vector& protein); + static bool isNewModel(const std::string & line); + + void clear(){_copies.clear();}; - size_t size() const {return _copies.size();}; +}; - AttractRigidbody& operator[](uint i){return _copies[i];}; -}; +class AttractMcop: public Mcop +{ +private: + + std::vector attract_copies; + +public: + + //Constructor: no arguments + AttractMcop(){ }; + //Constructor: takes filename + AttractMcop(std::string filename); + //Constructor: takes mcop + AttractMcop(const Mcop& mcop); + + virtual void addCopy(const AttractRigidbody& cop){attract_copies.push_back(cop);}; + + virtual size_t size() const {return attract_copies.size();}; + + virtual AttractRigidbody& operator[](uint i){return attract_copies[i];}; + virtual AttractRigidbody& getCopy(int i){return attract_copies[i];}; + +}; @@ -58,32 +81,47 @@ class Mcoprigid //multicopy rigidbody { public: - + //Constructor: no arguments Mcoprigid(); - + //Constructor: takes filename + Mcoprigid(std::string filename); //using default copy operator - void setMain(AttractRigidbody& main) ; - void addEnsemble(const Region& reg){ _vregion.push_back(reg); std::vector v; _weights.push_back(v); }; + void iniWeights(); + void setCore(AttractRigidbody& core) ; + void addEnsemble(const AttractMcop& reg){ _vregion.push_back(reg); std::vector v; _weights.push_back(v); }; void AttractEulerRotate(const dbl& phi, const dbl& ssi, const dbl& rot); void Translate(const Coord3D& c); - + Coord3D FindCenter(){return _core.FindCenter();}; void PrintWeights(); std::vector > getWeights(){return _weights;}; + + void ReadMcoprigidPDB(const std::string name); + void ReadMcoprigidPDB(std::istream& file, AttractRigidbody& core, std::vector& regions); + uint line_to_region_number(std::string line); + uint line_to_copy_number(std::string line); + + AttractRigidbody& getCore(){return _core;}; + std::vector& getRegions(){return _vregion;}; + AttractMcop& getRegion(int i){return _vregion[i];}; + size_t size() const {return _vregion.size();}; + virtual void setRotation(bool value) {_core.setRotation(value);}; ///< allow/disallow rotation + virtual void setTranslation(bool value) {_core.setTranslation(value);}; ///< allow/disallow translation + bool checkRotation(){return _core.checkRotation();}; + bool checkTranslation(){return _core.checkTranslation();}; private: - AttractRigidbody _main; - std::vector< Region > _vregion ; + AttractRigidbody _core; + std::vector< AttractMcop > _vregion ; - bool _complete ; - Coord3D _center ; ///
> _weights; + bool _complete ; //is this used ? + Coord3D _center ; ///
> _weights; friend class McopForceField; @@ -102,20 +140,24 @@ class McopForceField: public ForceField public: McopForceField(BaseAttractForceField& ff, dbl cutoff) - :_ff(ff), _cutoff(cutoff) {}; - + :_ff(ff), _cutoff(cutoff) { _beta = 0.3; }; + void ini_energies(); dbl Function(const Vdouble&); + dbl CalcEnergy(Mcoprigid & receptor,Mcoprigid & lig, BaseAttractForceField & ff, dbl cutoff); void Derivatives(const Vdouble& v, Vdouble & g ); + void setReceptor(Mcoprigid& rec); + void setLigand(Mcoprigid& lig); - void setReceptor(const Mcoprigid& rec) {_receptor = rec;}; - void setLigand(const Mcoprigid& lig) { _centered_ligand = lig; }; + uint ProblemSize(); + void initMinimization(){}; + virtual void saveWeights(){_savedWeights.push_back(_receptor.getWeights());}; + std::vector< std::vector< std::vector > > getSavedWeights(){return _savedWeights;}; - void calculate_weights(Mcoprigid& lig, bool print=false); + std::vector > getWeights(){return _receptor.getWeights();}; + std::vector > getMcopE(){return _Eik;}; - uint ProblemSize() {return 6;}; - void initMinimization(){}; private: @@ -127,6 +169,14 @@ class McopForceField: public ForceField Mcoprigid _centered_ligand ; Mcoprigid _moved_ligand ; Mcoprigid _receptor; + std::vector< std::vector > _Eik; //Loop copy interaction energies + std::vector< std::vector< std::vector > > _dEik; //- force applied to the ligand by each loop copy + dbl _E; //Total interaction energy + std::vector _Eregion; //Loop region interaction energies + std::vector _Z; // one Z per loop region + std::vector _Zprime; // one Z' per loop region + dbl _beta; // constant in the energy function + std::vector< std::vector< std::vector > > _savedWeights; //conserved weights throughout minimizations of a docking @@ -138,5 +188,3 @@ class McopForceField: public ForceField #endif // _MCOPFF_H_ - - diff --git a/headers/pdbio.h b/headers/pdbio.h index 7d76d00..43cd742 100644 --- a/headers/pdbio.h +++ b/headers/pdbio.h @@ -9,7 +9,6 @@ #include #include - #include "rigidbody.h" namespace PTools @@ -21,6 +20,9 @@ void ReadPDB(std::istream& fichier,Rigidbody& protein ); ///< read a PDB file fr void ReadPDB(const std::string name,Rigidbody& protein ); ///< read a PDB file from a filename and load data in Rigidbody void WritePDB(const Rigidbody& rigid, std::string filename); ///< write a PDB file given a Rigidbody and a filename +Coord3D pdbToCoords(const std::string & line); ///< returns Coord3D from a PDB-like string +void pdbToAtomproperty(const std::string & line, Atomproperty& atp); ///< Set Atomproperty from a PDB-like string +bool isAtom(const std::string & line ); ///< check if pdb line contains an atom } #endif //#ifndef PDBIO_H diff --git a/headers/rigidbody.h b/headers/rigidbody.h index c1ae3a4..38dc7ca 100644 --- a/headers/rigidbody.h +++ b/headers/rigidbody.h @@ -29,11 +29,10 @@ #include #include - +#include "coordsarray.h" #include "coord3d.h" #include "atom.h" #include "basetypes.h" -#include "coordsarray.h" namespace PTools @@ -44,7 +43,7 @@ namespace PTools class AtomSelection; // forward declaration -class Rigidbody:private CoordsArray +class Rigidbody: public CoordsArray { private: @@ -102,10 +101,10 @@ class Rigidbody:private CoordsArray void SetAtom(uint pos, const Atom& atom); /// add an atom to the molecule (deep copy) - void AddAtom(const Atomproperty& at, Coord3D co); + virtual void AddAtom(const Atomproperty& at, Coord3D co); /// add an atom to the molecule - void AddAtom(const Atom& at); + virtual void AddAtom(const Atom& at); //returns the coordinates of atom i Coord3D GetCoords(uint i) const diff --git a/src/attractforcefield.cpp b/src/attractforcefield.cpp index 2970e86..f9413d0 100644 --- a/src/attractforcefield.cpp +++ b/src/attractforcefield.cpp @@ -355,10 +355,11 @@ dbl BaseAttractForceField::Function(const Vdouble& stateVars ) //put the ligands to the correct positions defined by stateVars for (uint i=0; im_activeAtoms.size(); i++) { @@ -664,8 +662,11 @@ void BaseAttractForceField::Rota(uint molIndex, dbl phi,dbl ssi, dbl rot, Vdoubl Y = coords.y; Z = coords.z; + xar=X*crot+Y*srot; yar=-X*srot+Y*crot; + + pm[0][0]=-xar*cssp-yar*cp-Z*sssp ; pm[1][0]=xar*cscp-yar*sp+Z*sscp ; pm[2][0]=0.0 ; @@ -678,6 +679,7 @@ void BaseAttractForceField::Rota(uint molIndex, dbl phi,dbl ssi, dbl rot, Vdoubl pm[1][2]=yar*cssp-xar*cp ; pm[2][2]=-yar*ss ; + for (uint j=0;j<3;j++) { delta[j+shift] += pm[0][j] * pLigMoved->m_forces[atomIndex].x ; @@ -686,6 +688,7 @@ void BaseAttractForceField::Rota(uint molIndex, dbl phi,dbl ssi, dbl rot, Vdoubl } } + if (print) std::cout << "Rotational forces: " << delta[shift] << " " << delta[shift+1] << " " << delta[shift+2] << std::endl; return; diff --git a/src/attractrigidbody.cpp b/src/attractrigidbody.cpp index aa25d56..f130214 100644 --- a/src/attractrigidbody.cpp +++ b/src/attractrigidbody.cpp @@ -84,7 +84,27 @@ std::swap(m_activeAtoms, newactivelist); } +void AttractRigidbody::AddAtom(const Atom& at) +{ + Atomproperty atp(at); + Coord3D co = at.coords; + AddAtom(atp,co); +} +void AttractRigidbody::AddAtom(const Atomproperty& at, Coord3D co) +{ + mAtomProp.push_back(at); + AddCoord(co); + //Add AttractRigidbody parameters + uint atcategory = 0; + dbl atcharge = 0.0; + std::string extra = at.extra; + std::istringstream iss( extra ); + iss >> atcategory >> atcharge ; + m_atomTypeNumber.push_back(atcategory-1); // -1 to directly use the atomTypeNumber into C-array + m_charge.push_back(atcharge); + m_forces.push_back(Coord3D()); +} } //namespace PTools diff --git a/src/mcopff.cpp b/src/mcopff.cpp index 56b623c..a0ff212 100644 --- a/src/mcopff.cpp +++ b/src/mcopff.cpp @@ -1,23 +1,42 @@ #include "mcopff.h" #include - +#include +#include +#include namespace PTools { +/////////////////// -- Class Mcoprigid -- //////////////////// -Mcoprigid::Mcoprigid() -{ +Mcoprigid::Mcoprigid(){ _complete = false; }; +Mcoprigid::Mcoprigid(std::string filename){ + ReadMcoprigidPDB(filename); + _center = _core.FindCenter(); + iniWeights(); + +} +// Normalized weight initialized at 1/(number of copies) and denormalized weights initialized at 0. +void Mcoprigid::iniWeights(){ + for(uint i=0; i < _vregion.size(); i++){ + std::vector newvector; + _weights.push_back(newvector); + for(uint j=0; j < _vregion[i].size(); j++){ + dbl weight = 1.0/(double)_vregion[i].size(); + _weights[i].push_back(weight); + } + } +} -void Mcoprigid::setMain(AttractRigidbody& main) { - _main=main; - _center = _main.FindCenter(); +void Mcoprigid::setCore(AttractRigidbody& core) { + _core=core; + _center = _core.FindCenter(); }; @@ -27,10 +46,8 @@ void Mcoprigid::AttractEulerRotate(const dbl& phi, const dbl& ssi, const dbl& ro //Warning: makes euler rotation without centering //the Mcoprigid object must be centered - -//rotates the main body: - _main.AttractEulerRotate(phi, ssi, rot); - +//rotates the core body: + _core.AttractEulerRotate(phi, ssi, rot); //for each multicopy region, rotates the copy: for (uint i=0; i < _vregion.size(); i++) @@ -42,17 +59,73 @@ void Mcoprigid::AttractEulerRotate(const dbl& phi, const dbl& ssi, const dbl& ro void Mcoprigid::Translate(const Coord3D& c) { -//translates the main body: - _main.Translate(c); + +//translates the core body: + _core.Translate(c); //for each multicopy region, translates the copy: for (uint i=0; i < _vregion.size(); i++) - for (uint j=0; j<_vregion[i].size(); j++) + for (uint j=0; j < _vregion[i].size(); j++) _vregion[i][j].Translate(c); } +void Mcoprigid::ReadMcoprigidPDB(const std::string name) { + + // pointer toward the filename given in the constructor argument + std::ifstream file(name.c_str()); + if (!file) + { + std::ostringstream oss; + throw std::invalid_argument("##### ReadPDB:Could not open file \"" + name + "\" #####") ; + } + + ReadMcoprigidPDB(file, _core, _vregion); + file.close(); + +} + + +void Mcoprigid::ReadMcoprigidPDB(std::istream& file, AttractRigidbody& core, std::vector& regions){ + + uint region_num = 0; + uint copy_num = 0; + std::string line; + while(std::getline(file, line)){ + if(Mcop::isNewModel(line)){ + // The line is a new model + if(line.size() > 12) region_num = line_to_region_number(line); + if(line.size() > 15) copy_num = line_to_copy_number(line); + AttractRigidbody model; + while(std::getline(file,line)){ + if(isAtom(line)){ + // The line is an atom + Coord3D pos = pdbToCoords(line); + Atomproperty a; + pdbToAtomproperty(line, a); + // if region_num is core, append attom to core + if(region_num == 0) core.AddAtom(a,pos); + else model.AddAtom(a,pos); + } + + else{ + // if region_num is not core, add region copy + if(region_num != 0){ + // if new region, add region to vector before adding region copy + if(regions.size() < region_num) + regions.push_back(AttractMcop()); + // copy + regions[region_num-1].addCopy(model); + } + //Just finished adding a region copy + break; + } + } + } + } +} + void Mcoprigid::PrintWeights() { @@ -67,73 +140,161 @@ void Mcoprigid::PrintWeights() } } +uint Mcoprigid::line_to_region_number(std::string line){ + return std::atoi(line.substr(12,1).c_str()); +} +uint Mcoprigid::line_to_copy_number(std::string line){ + return std::atoi(line.substr(15,line.size()-15).c_str()); +} -/////////////////////////////////////////////////// -// Forcefield implementation -/////////////////////////////////////////////////// +/////////////////// -- Class Mcop -- //////////////////// -void McopForceField::calculate_weights(Mcoprigid& lig, bool print) -{ +Mcop::Mcop(std::string filename){ + ReadmcopPDB(filename); +} -//loop over copies regions - for (uint loopregion=0; loopregion < _receptor._vregion.size() ; loopregion++) +void Mcop::ReadmcopPDB(const std::string name) { + + std::string nomfich = name ; + // pointer toward the filename given in the constructor argument + std::ifstream file(nomfich.c_str()); + if (!file) { - if (print) - { - std::cout << " Region: " << loopregion << "\n"; - } + std::ostringstream oss; + throw std::invalid_argument("##### ReadPDB:Could not open file \"" + nomfich + "\" #####") ; + } - //calculates interaction energy between receptor copies and ligand body: - std::vector Eik; + ReadmcopPDB(file, _copies); + file.close(); - for (uint copy = 0; copy < _receptor._vregion[loopregion].size(); copy++) - { - AttractPairList cpl ( lig._main, _receptor._vregion[loopregion][copy], _cutoff ); - dbl e = _ff.nonbon8( lig._main, _receptor._vregion[loopregion][copy] , cpl ); - Eik.push_back(e); +} + +void Mcop::ReadmcopPDB(std::istream& file, std::vector& protein) { + + ReadModelsPDB(file, protein); + +} + +void Mcop::ReadModelsPDB(std::istream& file, std::vector& protein){ + + std::string line; + while(std::getline(file, line)){ + if(isNewModel(line)){ + // The line is a new model + Rigidbody model; + while(std::getline(file,line)){ + if(isAtom(line)){ + // The line is an atom + Coord3D pos = pdbToCoords(line); + Atomproperty a; + pdbToAtomproperty(line, a); + model.AddAtom(a,pos); + } + else{ + protein.push_back(model); + //Just finished adding a model + break; + } + } } + } +} +bool Mcop::isNewModel(const std::string & line){ - const dbl RT = 0.592 * 298.0 ; + if(line.substr(0,5)==(std::string)"MODEL") return true; + else return false; +} - //calculate weights - dbl sumweights = 0.0; - std::vector weights; - for (uint i=0; i newvector; + std::vector< std::vector > newcoordsvecvec; + std::vector threezeros(3, 0.0); + _Eik.push_back(newvector); + _Z.push_back(0); + _Zprime.push_back(0); + _Eregion.push_back(0); + _dEik.push_back(newcoordsvecvec); + for(uint k=0; k < _receptor.getRegions()[i].size(); k++){ + _Eik[i].push_back(0); + std::vector newcoordsvec; + _dEik[i].push_back(newcoordsvec); + } + } } +uint McopForceField::ProblemSize() +{ + + uint size = 0; + if (_centered_ligand.getCore().hastranslation) size += 3; + if (_centered_ligand.getCore().hasrotation) size += 3; + return size; +} @@ -144,133 +305,371 @@ void McopForceField::calculate_weights(Mcoprigid& lig, bool print) */ dbl McopForceField::Function(const Vdouble & v) { - - dbl ener = 0.0 ; - dbl enercopy =0.0; + dbl ener_region = 0.0 ; + dbl ener_core = 0.0 ; // 1) put the objects to the right place - + _moved_ligand = _centered_ligand; Mcoprigid & lig = _moved_ligand ; assert(lig._vregion.size()==0); - + //TODO: take into account if no rotation or no translation lig.AttractEulerRotate(v[0],v[1],v[2]); lig.Translate(Coord3D(v[3],v[4],v[5])); - + assert(_receptor._vregion.size() == _receptor._weights.size()); //2) calculates the energy - //2.1) main ligand body with main receptor + //2.1) core ligand body with core receptor + AttractPairList pl (_receptor._core, lig._core, _cutoff ); + ener_core += _ff.nonbon8(_receptor._core, lig._core, pl ); - AttractPairList pl (_receptor._main, lig._main, _cutoff ); - ener += _ff.nonbon8(_receptor._main, lig._main, pl ); + //2.2) core lignd with receptor copies: + std::vector< std::vector >& weights = _receptor._weights; + std::vector< std::vector >& Eik = _Eik; + std::vector< std::vector< std::vector > >& dEik = _dEik; + std::vector& Z = _Z; + std::vector& Zprime = _Zprime; + dbl& E = _E; + std::vector& Eregion = _Eregion; + dbl& beta = _beta; - //2.2) main lignd with receptor copies: - assert(_receptor._vregion.size() == _receptor._weights.size()); - for (uint loopregion=0; loopregion < _receptor._vregion.size() ; loopregion++) { //calculates interaction energy between receptor copies and ligand body: -// std::vector Eik; - + dbl enercopy =0.0; + dbl& Z_i = Z[loopregion]; + dbl& Zprime_i = Zprime[loopregion]; - Region& ref_ensemble = _receptor._vregion[loopregion]; - std::vector& ref_weights = _receptor._weights[loopregion]; + Z_i = 0.0; + Zprime_i = 0.0; - assert( ref_ensemble.size() == ref_weights.size()); + AttractMcop& ref_ensemble = _receptor._vregion[loopregion]; + std::vector& weights_loop = weights[loopregion]; + + assert( ref_ensemble.size() == weights_loop.size()); + double Emin = DBL_MAX; for (uint copynb = 0; copynb < ref_ensemble.size(); copynb++) { - dbl& weight = ref_weights[copynb]; + AttractRigidbody& copy = ref_ensemble[copynb]; - AttractPairList cpl ( lig._main, copy, _cutoff ); + AttractPairList cpl (copy, lig._core, _cutoff ); std::vector copyforce(copy.Size()); - std::vector mainforce(lig._main.Size()); + std::vector coreforce(lig._core.Size()); + + enercopy = _ff.nonbon8_forces(copy, lig._core, cpl, copyforce, coreforce); + Eik[loopregion][copynb] = enercopy; + if(enercopy < Emin) Emin = enercopy; + + //add force to ligand + assert(lig._core.Size() == coreforce.size()); + + dEik[loopregion][copynb] = coreforce; + } -// dbl e = _ff.nonbon8( lig._main, _receptor._vregion[loopregion][copy] , cpl ); - dbl e = _ff.nonbon8_forces(lig._main, copy, cpl, mainforce, copyforce); - enercopy += e * weight;//lig._weights[loopregion][copy]; + for (uint copynb = 0; copynb < ref_ensemble.size(); copynb++){ + Z_i += exp(-beta*(Eik[loopregion][copynb] - Emin)); + Zprime_i += -Eik[loopregion][copynb]*exp(-beta*(Eik[loopregion][copynb] - Emin)); + } - //multiply forces by copy weight: - for(uint i=0; i >& weights = receptor._weights; + std::vector< std::vector >& Eik = _Eik; + std::vector& Z = _Z; + std::vector& Zprime = _Zprime; + dbl& E = _E; + dbl beta = 1; + for (uint loopregion=0; loopregion < receptor._vregion.size() ; loopregion++) + { -void McopForceField::Derivatives(const Vdouble& v, Vdouble & g ) -{ + //calculates interaction energy between receptor copies and ligand body: + + dbl enercopy =0.0; + dbl& Z_i = Z[loopregion]; + dbl& Zprime_i = Zprime[loopregion]; -//sum the forces over x, y and z: + Z_i = 0.0; + Zprime_i = 0.0; -Coord3D ligtransForces; //translational forces for the ligand: -for(uint i=0; i<_moved_ligand._main.Size(); i++) - { - ligtransForces += _moved_ligand._main.m_forces[i]; - } + + AttractMcop& ref_ensemble = receptor._vregion[loopregion]; + std::vector& weights_loop = weights[loopregion]; + + assert( ref_ensemble.size() == weights_loop.size()); -Coord3D receptortransForces; -for(uint i=0; i<_receptor._main.Size(); i++) + for (uint copynb = 0; copynb < ref_ensemble.size(); copynb++) + { + + AttractRigidbody& copy = ref_ensemble[copynb]; + + AttractPairList cpl (copy, lig._core, cutoff ); + std::vector copyforce(copy.Size()); + std::vector coreforce(lig._core.Size()); + + enercopy = ff.nonbon8_forces(copy, lig._core, cpl, copyforce, coreforce); + Eik[loopregion][copynb] = enercopy; + + Z_i += exp(-beta*enercopy); + Zprime_i += -enercopy*exp(-beta*enercopy); + + } + for (uint copynb = 0; copynb < ref_ensemble.size(); copynb++) + { + dbl& weight = weights_loop[copynb]; + weight = exp(-beta*Eik[loopregion][copynb])/Z_i; + } + ener_region += -Zprime_i/Z_i; + } + E = ener_core + ener_region; + return E; +} + + + +void McopForceField::Derivatives(const Vdouble& v, Vdouble & g ) { -receptortransForces+= _receptor._main.m_forces[i]; + +Mcoprigid & lig = _moved_ligand; +uint svptr = 0; // stateVars 'pointer' +// - force applied to the ligand by each loop copy +std::vector< std::vector< std::vector > >& dEik = _dEik; +//Interaction nergy between each loop copy and the ligand +std::vector< std::vector >& Eik = _Eik; +//Interactino energy between each variable region (sum of its copies) and the ligand. +std::vector& Eregion = _Eregion; +//Weight of each loop copy +std::vector< std::vector > & weights = _receptor._weights; +//beta coefficient (constant in energy function) +dbl& beta = _beta; + +std::vector< std::vector > coef; +std::vector newvector; +// iterate through loop regions +for(uint i=0; i < dEik.size(); i++){ + coef.push_back(newvector); + // iterate through loop copies + for(uint k=0; k < dEik[i].size(); k++){ + coef[i].push_back((1+beta*(Eregion[i] - Eik[i][k]))*weights[i][k]); + } } +// calculate de rotational forces: +if (lig.getCore().hasrotation){ + //printf("hasrotation\n"); + dbl cs,cp,ss,sp,cscp,sscp,sssp,crot,srot,xar,yar,cssp,x,y,z; + dbl pm[3][3]; -for (uint i=0; i <_receptor._vregion.size(); i++) - { - Region& ens = _receptor._vregion[i]; - std::vector & weights = _receptor._weights[i]; + // !c calculates orientational force contributions + // !c component 1: phi-angle + // !c component 2: ssi-angle + // !c component 3: rot-angle - for(uint j=0; j dEcore(3, 0.0); + std::vector dEloops(3, 0.0); + + for(uint atomIndex=0; atomIndex<_moved_ligand._core.Size(); atomIndex++){ + + Coord3D coords = pLigCentered->GetCoords(atomIndex); + + x = coords.x; + y = coords.y; + z = coords.z; + + xar=x*crot+y*srot; + yar=-x*srot+y*crot; + + pm[0][0]=-xar*cssp-yar*cp-z*sssp ; + pm[1][0]=xar*cscp-yar*sp+z*sscp ; + pm[2][0]=0.0 ; + + pm[0][1]=-xar*sscp+z*cscp ; + pm[1][1]=-xar*sssp+z*cssp ; + pm[2][1]=-xar*cs-z*ss ; + + pm[0][2]=yar*cscp+xar*sp ; + pm[1][2]=yar*cssp-xar*cp ; + pm[2][2]=-yar*ss ; + + //adding the derivatives for the core: + for(uint j=0; j < 3; j++){ + dEcore[j] += pm[0][j] * pLigMoved->m_forces[atomIndex].x; + dEcore[j] += pm[1][j] * pLigMoved->m_forces[atomIndex].y; + dEcore[j] += pm[2][j] * pLigMoved->m_forces[atomIndex].z; + } - } + // iterate through loop regions + for(uint i=0; i < dEik.size(); i++){ + // iterate through loop copies + for(uint k=0; k < dEik[i].size(); k++){ + // iterate through the three rotation angles + for(uint j=0; j < 3; j++){ + dEloops[j] += coef[i][k] * pm[0][j] * dEik[i][k][atomIndex].x; + dEloops[j] += coef[i][k] * pm[1][j] * dEik[i][k][atomIndex].y; + dEloops[j] += coef[i][k] * pm[2][j] * dEik[i][k][atomIndex].z; + } + } + } + } + for(uint j=0; j < 3; j++){ + g[svptr+j] += dEcore[j] + dEloops[j]; + } + svptr += 3; } +//sum the forces over x, y and z: calculate de translational forces +if (lig.getCore().hastranslation){ + + g[svptr+0] = 0.0; + g[svptr+1] = 0.0; + g[svptr+2] = 0.0; + + dbl flim = 1.0e18; + dbl ftr1, ftr2, ftr3, fbetr; + ftr1=0.0; + ftr2=0.0; + ftr3=0.0; + + Coord3D ligtransForces; //translational forces for the ligand from loopregion + + //initializing derivative variables + std::vector dEcore(3, 0.0); + std::vector dEloops(3, 0.0); + + //adding derivatives + for(uint atomIndex=0; atomIndex<_moved_ligand._core.Size(); atomIndex++){ + ligtransForces += _moved_ligand._core.m_forces[atomIndex]; + for(uint i=0; i < dEik.size(); i++){ + for(uint k=0; k < dEik[i].size(); k++){ + dEloops[0] += coef[i][k] * dEik[i][k][atomIndex].x; + dEloops[1] += coef[i][k] * dEik[i][k][atomIndex].y; + dEloops[2] += coef[i][k] * dEik[i][k][atomIndex].z; + } + } + } + dEcore[0] = ligtransForces.x; + dEcore[1] = ligtransForces.y; + dEcore[2] = ligtransForces.z; + ftr1 = dEcore[0] + dEloops[0]; + ftr2 = dEcore[1] + dEloops[1]; + ftr3 = dEcore[2] + dEloops[2]; + //ftr1 =g[svptr+0]; + //ftr2 =g[svptr+1]; + //ftr3 =g[svptr+2]; + + // force reduction, some times helps in case of very "bad" start structure + for (uint i=0; i<3; i++) + { + fbetr=ftr1*ftr1 +ftr2*ftr2 +ftr3*ftr3; + if (fbetr > flim) + { + ftr1=.01*ftr1; + ftr2=.01*ftr2; + ftr3=.01*ftr3; + } + } + + + assert(svptr+2 < g.size()); + g[svptr+0] = ftr1; + g[svptr+1] = ftr2; + g[svptr+2] = ftr3; + svptr += 3; +} + + +} + } // namespace PTools + + + + + diff --git a/src/minimizers/lbfgs_interface.cpp b/src/minimizers/lbfgs_interface.cpp index a658d95..c592740 100644 --- a/src/minimizers/lbfgs_interface.cpp +++ b/src/minimizers/lbfgs_interface.cpp @@ -1,4 +1,5 @@ #include "lbfgs_interface.h" +#include "mcopff.h" #include #include @@ -71,18 +72,15 @@ inline std::vector todbl(std::vector & vcplx) void Lbfgs::minimize(int maxiter) { - int n = objToMinimize.ProblemSize(); std::cout << "number of free variables for the minimizer: " << n << std::endl; - - std::vector l(n); std::vector u(n); Vint nbd(n); - x.resize(n); g.resize(n); + for (int i=0;iiprint=-1; double f = DBL_MAX; @@ -109,8 +105,10 @@ void Lbfgs::minimize(int maxiter) m_opt->max_iter = maxiter; int last_iter = 0; + int sub_iter = 0; + + double grad = 0; - /* opt->iprint = 0;*/ while (1) { rc = lbfgsb_run(m_opt, &x[0], &f, &g[0]); if (rc == 0) { @@ -120,65 +118,31 @@ void Lbfgs::minimize(int maxiter) break; } else if (rc == 1) { - - -/* -//check analytical derivatives with surreal: -{ - std::vector vdblx; - tocplx(x,vdblx); - std::vector vdblg; - tocplx(g,vdblg); - - f = objToMinimize.Function(vdblx); - objToMinimize.Derivatives(vdblx,vdblg); - - g=todbl(vdblg); - - - for (uint i=0; i svdblx = vdblx ; - svdblx[i]+=surreal(0,1); - std::cout <<"check: " << svdblx[i] << std::endl; - - std::cout << g[i] << "==" << imag(objToMinimize.Function(svdblx)) << " " ; - } - std::cout << std::endl; -} -//end check derivatives -*/ - - - - if (last_iter < m_opt->niter) { //this is a new iteration last_iter = m_opt->niter; + sub_iter = 0; //saves the minimizer variables for each iteration (can be useful for generating animations) m_vars_over_time.push_back(x); - + objToMinimize.saveWeights(); //Only saves weights if objToMinimize is McopForceField. } - + + double sum = 0; + for(int i=0; i < g.size(); i++){ + sum += g[i]*g[i]; + } + std::vector vdblx; tocplx(x,vdblx); std::vector vdblg; tocplx(g,vdblg); - f = objToMinimize.Function(vdblx); objToMinimize.Derivatives(vdblx,vdblg); - g=todbl(vdblg); - -// std::cout << "analytical derivatives: \n"; -// for(uint i=0; i=m_vars_over_time.size()) return m_vars_over_time[iter]; } - - +std::vector< std::vector > Lbfgs::getWeights(){ + if (McopForceField * p = dynamic_cast(&objToMinimize)){ + // objToMinimize is of type McopForceField + ForceField& r_objToMinimize = objToMinimize; + McopForceField& r_Mcop_objToMinimize = dynamic_cast(r_objToMinimize); + r_Mcop_objToMinimize.getWeights(); + } +} diff --git a/src/rigidbody.cpp b/src/rigidbody.cpp index cfa71c0..2c6510b 100644 --- a/src/rigidbody.cpp +++ b/src/rigidbody.cpp @@ -144,7 +144,9 @@ dbl Rigidbody::Radius() void Rigidbody::Translate(const Coord3D& tr) { + CoordsArray::Translate(tr); + } void Rigidbody::AttractEulerRotate(dbl phi, dbl ssi, dbl rot)