Python molecule.Molecule类代码示例

 def read_quantum():
     Result = None
     os.chdir('calcs')
     for i in range(M.ns):
         dnm = eval(formstr % i)
         print("\rNow in directory %i" % i, end=' ')
         if os.path.exists(dnm):
             os.chdir(dnm)
             if os.path.exists('qchem.out'):
                 Output = Molecule('qchem.out')
                 if os.path.exists('plot.esp'):
                     ESP = Molecule('plot.esp')
                     #print ESP.Data.keys()
                     Output.qm_espxyzs = list(ESP.qm_espxyzs)
                     Output.qm_espvals = list(ESP.qm_espvals)
                     #Output += Molecule('plot.esp')
                 if Result == None:
                     Result = Output
                 else:
                     Result += Output
             else:
                 raise Exception("The output file %s doesn't exist." % os.path.abspath('qchem.out'))
             os.chdir('..')
         else:
             raise Exception("The subdirectory %s doesn't exist." % os.path.abspath(dnm))
     os.chdir('..')
     return Result

def QChem_Dielectric_Energy(fnm,wq):
    QCIn = Molecule(fnm)
    for i in range(QCIn.na):
        # Q-Chem crashes if it doesn't recognize the chemical element
        if QCIn.Data['elem'][i] in ['M','L']:
            QCIn.Data['elem'][i] = 'He'
    CalcDir=os.path.splitext(fnm)[0]+".d"
    GoInto(CalcDir)
    digits = len(str(QCIn.ns))
    for i in range(QCIn.ns):
        sdir = "%%0%ii" % digits % i
        GoInto(sdir)
        QCIn.write("qchem.in",select=i)
        queue_up(wq,"qchem40 qchem.in qchem.out",input_files=["qchem.in"],output_files=["qchem.out"],verbose=False)
        Leave(sdir)
    wq_wait(wq,verbose=False)
    PCM_Energies = []
    for i in range(QCIn.ns):
        sdir = "%%0%ii" % digits % i
        GoInto(sdir)
        for line in open("qchem.out"):
            if "PCM electrostatic energy" in line:
                PCM_Energies.append(float(line.split()[-2]))
        Leave(sdir)
    Leave(CalcDir)
    return np.array(PCM_Energies) * 2625.5

def main():
    M = Molecule(argv[1])
    tempfnm = argv[3] if len(argv) >= 4 else None    
    if tempfnm != None:
        M.add_quantum(tempfnm)
    print M.Data.keys()
    M.write(argv[2])

    def __init__(self,options,tgt_opts,forcefield):
        super(BindingEnergy,self).__init__(options,tgt_opts,forcefield)
        self.set_option(None, None, 'inter_txt', os.path.join(self.tgtdir,tgt_opts['inter_txt']))
        self.global_opts, self.sys_opts, self.inter_opts = parse_interactions(self.inter_txt)
        # If the global option doesn't exist in the system / interaction, then it is copied over.
        for opt in self.global_opts:
            for sys in self.sys_opts:
                if opt not in self.sys_opts[sys]:
                    self.sys_opts[sys][opt] = self.global_opts[opt]
            for inter in self.inter_opts:
                if opt not in self.inter_opts[inter]:
                    self.inter_opts[inter][opt] = self.global_opts[opt]
        for inter in self.inter_opts:
            if 'energy_unit' in self.inter_opts[inter] and self.inter_opts[inter]['energy_unit'].lower() not in ['kilocalorie_per_mole', 'kilocalories_per_mole']:
                logger.error('Usage of physical units is has been removed, please provide all binding energies in kcal/mole\n')
                raise RuntimeError
            self.inter_opts[inter]['reference_physical'] = self.inter_opts[inter]['energy']

        if tgt_opts['energy_denom'] == 0.0:
            self.set_option(None, None, 'energy_denom', val=np.std(np.array([val['reference_physical'] for val in self.inter_opts.values()])))
        else:
            self.set_option(None, None, 'energy_denom', val=tgt_opts['energy_denom'])

        self.set_option(None, None, 'rmsd_denom', val=tgt_opts['rmsd_denom'])

        self.set_option(tgt_opts,'cauchy')
        self.set_option(tgt_opts,'attenuate')
        ## LPW 2018-02-11: This is set to True if the target calculates
        ## a single-point property over several existing snapshots.
        self.loop_over_snapshots = False

        logger.info("The energy denominator is: %s\n" % str(self.energy_denom)) 
        logger.info("The RMSD denominator is: %s\n" % str(self.rmsd_denom))

        if self.cauchy:
            logger.info("Each contribution to the interaction energy objective function will be scaled by 1.0 / ( denom**2 + reference**2 )\n")
            if self.attenuate:
                logger.error('attenuate and cauchy are mutually exclusive\n')
                raise RuntimeError
        elif self.attenuate:
            logger.info("Repulsive interactions beyond energy_denom will be scaled by 1.0 / ( denom**2 + (reference-denom)**2 )\n")
        ## Build keyword dictionaries to pass to engine.
        engine_args = OrderedDict(list(self.OptionDict.items()) + list(options.items()))
        del engine_args['name']
        ## Create engine objects.
        self.engines = OrderedDict()
        for sysname,sysopt in self.sys_opts.items():
            M = Molecule(os.path.join(self.root, self.tgtdir, sysopt['geometry']))
            if 'select' in sysopt:
                atomselect = np.array(uncommadash(sysopt['select']))
                M = M.atom_select(atomselect)
            if self.FF.rigid_water: M.rigid_water()
            self.engines[sysname] = self.engine_(target=self, mol=M, name=sysname, tinker_key=os.path.join(sysopt['keyfile']), **engine_args)

    def optimize(self, shot=0, method="newton", crit=1e-4):

        """ Optimize the geometry and align the optimized geometry to the starting geometry. """

        if os.path.exists("%s.xyz_2" % self.name):
            os.unlink("%s.xyz_2" % self.name)

        self.mol[shot].write("%s.xyz" % self.name, ftype="tinker")

        if method == "newton":
            if self.rigid:
                optprog = "optrigid"
            else:
                optprog = "optimize"
        elif method == "bfgs":
            if self.rigid:
                optprog = "minrigid"
            else:
                optprog = "minimize"

        o = self.calltinker("%s %s.xyz %f" % (optprog, self.name, crit))
        # Silently align the optimized geometry.
        M12 = Molecule("%s.xyz" % self.name, ftype="tinker") + Molecule("%s.xyz_2" % self.name, ftype="tinker")
        if not self.pbc:
            M12.align(center=False)
        M12[1].write("%s.xyz_2" % self.name, ftype="tinker")
        rmsd = M12.ref_rmsd(0)[1]
        cnvgd = 0
        mode = 0
        for line in o:
            s = line.split()
            if len(s) == 0:
                continue
            if "Optimally Conditioned Variable Metric Optimization" in line:
                mode = 1
            if "Limited Memory BFGS Quasi-Newton Optimization" in line:
                mode = 1
            if mode == 1 and isint(s[0]):
                mode = 2
            if mode == 2:
                if isint(s[0]):
                    E = float(s[1])
                else:
                    mode = 0
            if "Normal Termination" in line:
                cnvgd = 1
        if not cnvgd:
            for line in o:
                logger.info(str(line) + "\n")
            logger.info("The minimization did not converge in the geometry optimization - printout is above.\n")
        return E, rmsd

    def readsrc(self, **kwargs):

        """ Called by __init__ ; read files from the source directory. """

        self.leapcmd = onefile(kwargs.get('leapcmd'), 'leap', err=True)
        self.absleap = os.path.abspath(self.leapcmd)

        # Name of the molecule, currently just call it a default name.
        self.mname = 'molecule'

        if 'mol' in kwargs:
            self.mol = kwargs['mol']
        elif 'coords' in kwargs:
            crdfile = onefile(kwargs.get('coords'), None, err=True)
            self.mol = Molecule(crdfile, build_topology=False)

        # AMBER has certain PDB requirements, so we will absolutely require one.
        needpdb = True
        # if hasattr(self, 'mol') and all([i in self.mol.Data.keys() for i in ["chain", "atomname", "resid", "resname", "elem"]]):
        #     needpdb = False

        # Determine the PDB file name.
        # If 'pdb' is provided to Engine initialization, it will be used to 
        # copy over topology information (chain, atomname etc.).  If mol/coords
        # is not provided, then it will also provide the coordinates.
        pdbfnm = onefile(kwargs.get('pdb'), 'pdb' if needpdb else None, err=needpdb)
        if pdbfnm != None:
            mpdb = Molecule(pdbfnm, build_topology=False)
            if hasattr(self, 'mol'):
                for i in ["chain", "atomname", "resid", "resname", "elem"]:
                    self.mol.Data[i] = mpdb.Data[i]
            else:
                self.mol = copy.deepcopy(mpdb)
        self.abspdb = os.path.abspath(pdbfnm)

    def readsrc(self, **kwargs):

        """ Called by __init__ ; read files from the source directory. """

        self.key = onefile('key', kwargs['tinker_key'] if 'tinker_key' in kwargs else None)
        self.prm = onefile('prm', kwargs['tinker_prm'] if 'tinker_prm' in kwargs else None)
        if 'mol' in kwargs:
            self.mol = kwargs['mol']
        elif 'coords' in kwargs:
            if kwargs['coords'].endswith('.xyz'):
                self.mol = Molecule(kwargs['coords'], ftype="tinker")
            else:
                self.mol = Molecule(kwargs['coords'])
        else:
            arcfile = onefile('arc')
            if not arcfile: raise RuntimeError('Cannot determine which .arc file to use')
            self.mol = Molecule(arcfile)

    def readsrc(self, **kwargs):

        """ Called by __init__ ; read files from the source directory. """

        self.key = onefile("key", kwargs["tinker_key"] if "tinker_key" in kwargs else None)
        self.prm = onefile("prm", kwargs["tinker_prm"] if "tinker_prm" in kwargs else None)
        if "mol" in kwargs:
            self.mol = kwargs["mol"]
        elif "coords" in kwargs:
            if kwargs["coords"].endswith(".xyz"):
                self.mol = Molecule(kwargs["coords"], ftype="tinker")
            else:
                self.mol = Molecule(kwargs["coords"])
        else:
            arcfile = onefile("arc")
            if not arcfile:
                raise RuntimeError("Cannot determine which .arc file to use")
            self.mol = Molecule(arcfile)

def gather_generations():
    shots = Molecule('shots.gro')
    qdata = Molecule('qdata.txt')
    A1    = np.array(shots.xyzs)
    A2    = np.array(qdata.xyzs)
    if A1.shape != A2.shape:
        raise Exception('shots.gro and qdata.txt appear to contain different data')
    elif np.max(np.abs((A1 - A2).flatten())) > 1e-4:
        raise Exception('shots.gro and qdata.txt appear to contain different xyz coordinates')
    shots.qm_energies = qdata.qm_energies
    shots.qm_forces   = qdata.qm_forces
    shots.qm_espxyzs     = qdata.qm_espxyzs
    shots.qm_espvals     = qdata.qm_espvals
    First = True
    if First:
        All = shots
    else:
        All += shots
    return All

 def prepare_temp_directory(self, options, tgt_opts):
     abstempdir = os.path.join(self.root,self.tempdir)
     if self.FF.rigid_water:
         self.optprog = "optrigid"
         #LinkFile(os.path.join(self.root,self.tgtdir,"rigid.key"),os.path.join(abstempdir,"rigid.key"))
     else:
         self.optprog = "optimize"
     # Link the necessary programs into the temporary directory
     LinkFile(os.path.join(options['tinkerpath'],"analyze"),os.path.join(abstempdir,"analyze"))
     LinkFile(os.path.join(options['tinkerpath'],self.optprog),os.path.join(abstempdir,self.optprog))
     LinkFile(os.path.join(options['tinkerpath'],"superpose"),os.path.join(abstempdir,"superpose"))
     # Link the run parameter file
     # The master file might be unneeded??
     for sysname,sysopt in self.sys_opts.items():
         if self.FF.rigid_water:
             # Make every water molecule rigid
             # WARNING: Hard coded values here!
             M = Molecule(os.path.join(self.root, self.tgtdir, sysopt['geometry']),ftype="tinker")
             for a in range(0, len(M.xyzs[0]), 3):
                 flex = M.xyzs[0]
                 wat = flex[a:a+3]
                 com = wat.mean(0)
                 wat -= com
                 o  = wat[0]
                 h1 = wat[1]
                 h2 = wat[2]
                 r1 = h1 - o
                 r2 = h2 - o
                 r1 /= Np.linalg.norm(r1)
                 r2 /= Np.linalg.norm(r2)
                 # Obtain unit vectors.
                 ex = r1 + r2
                 ey = r1 - r2
                 ex /= Np.linalg.norm(ex)
                 ey /= Np.linalg.norm(ey)
                 Bond = 0.9572
                 Ang = Np.pi * 104.52 / 2 / 180
                 cosx = Np.cos(Ang)
                 cosy = Np.sin(Ang)
                 h1 = o + Bond*ex*cosx + Bond*ey*cosy
                 h2 = o + Bond*ex*cosx - Bond*ey*cosy
                 rig = Np.array([o, h1, h2]) + com
                 M.xyzs[0][a:a+3] = rig
             M.write(os.path.join(abstempdir,sysopt['geometry']),ftype="tinker")
         else:
             M = Molecule(os.path.join(self.root, self.tgtdir, sysopt['geometry']),ftype="tinker")
             if 'select' in sysopt:
                 atomselect = Np.array(uncommadash(sysopt['select']))
                 #atomselect = eval("Np.arange(M.na)"+sysopt['select'])
                 M = M.atom_select(atomselect)
             M.write(os.path.join(abstempdir,sysname+".xyz"),ftype="tinker")
             write_key_with_prm(os.path.join(self.root,self.tgtdir,sysopt['keyfile']),os.path.join(abstempdir,sysname+".key"),ffobj=self.FF)

    def readsrc(self, **kwargs):

        """ Called by __init__ ; read files from the source directory. """

        self.key = onefile(kwargs.get("tinker_key"), "key")
        self.prm = onefile(kwargs.get("tinker_prm"), "prm")
        if "mol" in kwargs:
            self.mol = kwargs["mol"]
        else:
            crdfile = onefile(kwargs.get("coords"), "arc", err=True)
            self.mol = Molecule(crdfile)

def main():
    M = Molecule(opts.coords)
    if len(M.molecules) != 1:
        raise RuntimeError('Input coordinates must be a single contiguous molecule')
    if opts.phi1 == None:
        raise RuntimeError('phi1 (the first quartet of atoms) must be provided')
    xyzout = []
    commout = []
    xyzcsh = []
    commcsh = []
    if opts.phi2 != None: # Two dimensional scan
        for inc1 in range(0, 360, opts.scan):
            for inc2 in range(0, 360, opts.scan):
                xyzrot, clash = get_rotated_xyz(M, [opts.phi1, opts.phi2], [inc1, inc2])
                print(inc1, inc2, "Clash" if clash else "Ok")
                comm = "Dihedrals %s, %s set to %i, %i" % (str(opts.phi1), str(opts.phi2), inc1, inc2)
                if clash: 
                    xyzcsh.append(xyzrot.copy())
                    commcsh.append(comm)
                else: 
                    xyzout.append(xyzrot.copy())
                    commout.append(comm)
    else: # One dimensional scan
        for inc1 in range(0, 360, opts.scan):
            xyzrot, clash = get_rotated_xyz(M, [opts.phi1], [inc1])
            print(inc1, "Clash" if clash else "Ok")
            comm = "Dihedral %s set to %i" % (str(opts.phi1), inc1)
            if clash: 
                xyzcsh.append(xyzrot.copy())
                commcsh.append(comm)
            else: 
                xyzout.append(xyzrot.copy())
                commout.append(comm)
    if len(xyzout) > 0:
        M.xyzs = xyzout
        M.comms = commout
        M.write(os.path.splitext(opts.coords)[0]+"_out"+os.path.splitext(opts.coords)[1])
    if len(xyzcsh) > 0:
        M.xyzs = xyzcsh
        M.comms = commcsh
        M.write(os.path.splitext(opts.coords)[0]+"_clash"+os.path.splitext(opts.coords)[1])

def main():
    # create molecule
    m = Molecule(opts.xyz_file)
    # write to pdb
    m.write(tf)
    # read into mdtraj
    p = mdtraj.load(tf)
    # get list of hydrogens to remove
    c_ndx = [i.index for i in p.topology.atoms if i.name in c_strip]
    tail_atms = [y for (x, y) in m.Data['bonds'] if x in c_ndx]
    tail_h = [i for i in tail_atms if i not in c_ndx]
    # throw out some hydrogens
    atoms_to_keep = [b.index for b in p.topology.atoms if b.index not in tail_h]
    p.restrict_atoms(atoms_to_keep)
    p.save(op)
    # hella circular 
    x = Molecule(op)
    x.write(opts.out_gro)

import numpy as np
from forcebalance.molecule import Molecule
from forcebalance.readfrq import read_frq_gen

# Frequency output file.
fout = sys.argv[1]

# Mode number, starting from 1.
modenum = int(sys.argv[2])

if modenum == 0:
    raise RuntimeError("Start mode number from one, please")

frqs, modes, intens, elem, xyz = read_frq_gen(fout)

M = Molecule()
M.elem = elem[:]
M.xyzs = []

xmode = modes[modenum - 1]
xmode /= (np.linalg.norm(xmode)/np.sqrt(M.na))
xmode *= 0.3 # Reasonable vibrational amplitude

spac = np.linspace(0, 1, 101)
disp = np.concatenate((spac, spac[::-1][1:], -1*spac[1:], -1*spac[::-1][1:-1]))

for i in disp:
    M.xyzs.append(xyz+i*xmode.reshape(-1,3))

M.comms = ['Vibrational Mode %i Frequency %.3f Displacement %.3f' % (modenum, frqs[modenum-1], disp[i]*(np.linalg.norm(xmode)/np.sqrt(M.na))) for i in range(len(M))]

#!/usr/bin/env python

from forcebalance.molecule import Molecule
from forcebalance.nifty import _exec
import os

np=""
if "CORES_PER_WORKER" in os.environ and int(os.environ["CORES_PER_WORKER"]) > 1:
    np=" -np %i" % int(os.environ["CORES_PER_WORKER"])

_exec("touch opt.xyz")
_exec("touch energy.txt")
_exec("rm -f qchem.out.prev")
_exec("touch qchem.out.prev")
qcin = Molecule("qchem.in", ftype="qcin")
qcin.edit_qcrems({'geom_opt_max_cycles':'100'})
qcin.write("qchem.in")
_exec("qchem42 %s qchem.in qchem.out &> qchem.err" % np)

def special_criterion():
    mk = 0
    mx = 0
    Cnvgd = {}
    for ln, line in enumerate(open("qchem.out").readlines()):
        if "Maximum optimization cycles reached" in line:
            mx = 1
        if "Maximum     Tolerance    Cnvgd?" in line:
            mk = ln
        if mk > 0 and ln > mk and ln <= mk+3:
            s = line.split()
            try: