当前位置: 首页>>代码示例>>Python>>正文


Python Molecule.write方法代码示例

本文整理汇总了Python中forcebalance.molecule.Molecule.write方法的典型用法代码示例。如果您正苦于以下问题:Python Molecule.write方法的具体用法?Python Molecule.write怎么用?Python Molecule.write使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在forcebalance.molecule.Molecule的用法示例。


在下文中一共展示了Molecule.write方法的14个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。

示例1: QChem_Dielectric_Energy

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
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
开发者ID:rmcgibbo,项目名称:forcebalance,代码行数:28,代码来源:qchemio.py

示例2: main

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
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])
开发者ID:JNapoli,项目名称:forcebalance,代码行数:9,代码来源:filecnv.py

示例3: prepare_temp_directory

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
 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)
开发者ID:rmcgibbo,项目名称:forcebalance,代码行数:54,代码来源:tinkerio.py

示例4: main

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
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)
开发者ID:kmckiern,项目名称:scripts,代码行数:20,代码来源:aa_to_crm_ua.py

示例5: main

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
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])
开发者ID:leeping,项目名称:forcebalance,代码行数:43,代码来源:dscan.py

示例6: int

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
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))]

M.write(os.path.splitext(fout)[0]+'.mode%03i.xyz' % modenum)
开发者ID:leeping,项目名称:forcebalance,代码行数:32,代码来源:anifrq.py

示例7: Molecule

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
#!/usr/bin/env python

from builtins import range
from forcebalance.molecule import Molecule
import os, sys

# Load in the Gromacs .gro file to be converted to TINKER format.
M = Molecule(sys.argv[1])

# Build the line suffix for the TINKER format.
tinkersuf = []
for i in range(M.na):
    if i%3==0:
        tinkersuf.append("%5i %5i %5i" % (1, i+2, i+3))
    else:
        tinkersuf.append("%5i %5i" % (2, i-i%3+1))
M.tinkersuf = tinkersuf

# Delete the periodic box.
del M.Data['boxes']

# Write the TINKER output format.
M.write(os.path.splitext(sys.argv[1])[0]+'.xyz', ftype='tinker')
开发者ID:leeping,项目名称:forcebalance,代码行数:25,代码来源:gro2arc-amoeba-water.py

示例8: Lipid

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]

#.........这里部分代码省略.........
                pt_label = "IC/%sK-%s%s" % (pt[0], pt[1], pt[2])
                if not os.path.exists(os.path.join(self.root, self.tgtdir, pt_label, self.lipid_coords)):
                    raise RuntimeError("Initial condition files don't exist; please provide IC directory")
                # Create molecule object for each IC.
                all_ic = Molecule(os.path.join(self.root, self.tgtdir, pt_label, self.lipid_coords))
                self.lipid_mols[pt] = []
                n_uniq_ic = int(self.RefData['n_ic'][pt])
                if n_uniq_ic > len(all_ic):
                    raise RuntimeError("Number of frames in initial conditions .gro file is less than the number of parallel simulations requested in data.csv")
                # Index ICs by pressure and temperature in a dictionary.
                for ic in range(n_uniq_ic):
                    self.lipid_mols[pt].append(all_ic[ic])
        else:
            # Read in lipid starting coordinates.
            if not os.path.exists(os.path.join(self.root, self.tgtdir, self.lipid_coords)): 
                logger.error("%s doesn't exist; please provide lipid_coords option\n" % self.lipid_coords)
                raise RuntimeError
            self.lipid_mol = Molecule(os.path.join(self.root, self.tgtdir, self.lipid_coords), toppbc=True)
            # Extra files to be linked into the temp-directory.
            self.nptfiles += [self.lipid_coords]
        # Scripts to be copied from the ForceBalance installation directory.
        self.scripts += ['npt_lipid.py']
        # Prepare the temporary directory.
        self.prepare_temp_directory()
        # Build keyword dictionary to pass to engine.
        if self.do_self_pol:
            self.gas_engine_args.update(self.OptionDict)
            self.gas_engine_args.update(options)
            del self.gas_engine_args['name']
            # Create engine object for gas molecule to do the polarization correction.
            self.gas_engine = self.engine_(target=self, mol=self.gas_mol, name="selfpol", **self.gas_engine_args)
        # Don't read indicate.log when calling meta_indicate()
        self.read_indicate = False
        self.write_indicate = False
        # Don't read objective.p when calling meta_get()
        self.read_objective = False
        #======================================#
        #          UNDER DEVELOPMENT           #
        #======================================#
        # Put stuff here that I'm not sure about. :)
        np.set_printoptions(precision=4, linewidth=100)
        np.seterr(under='ignore')
        ## Saved force field mvals for all iterations
        self.SavedMVal = {}
        ## Saved trajectories for all iterations and all temperatures
        self.SavedTraj = defaultdict(dict)
        ## Evaluated energies for all trajectories (i.e. all iterations and all temperatures), using all mvals
        self.MBarEnergy = defaultdict(lambda:defaultdict(dict))

    def prepare_temp_directory(self):
        """ Prepare the temporary directory by copying in important files. """
        abstempdir = os.path.join(self.root,self.tempdir)
        for f in self.nptfiles:
            LinkFile(os.path.join(self.root, self.tgtdir, f), os.path.join(abstempdir, f))
        for f in self.scripts:
            LinkFile(os.path.join(os.path.split(__file__)[0],"data",f),os.path.join(abstempdir,f))

    def read_data(self):
        # Read the 'data.csv' file. The file should contain guidelines.
        with open(os.path.join(self.tgtdir,'data.csv'),'rU') as f: R0 = list(csv.reader(f))
        # All comments are erased.
        R1 = [[sub('#.*$','',word) for word in line] for line in R0 if len(line[0]) > 0 and line[0][0] != "#"]
        # All empty lines are deleted and words are converted to lowercase.
        R = [[wrd.lower() for wrd in line] for line in R1 if any([len(wrd) for wrd in line]) > 0]
        global_opts = OrderedDict()
        found_headings = False
开发者ID:leeping,项目名称:forcebalance,代码行数:70,代码来源:lipid.py

示例9: AMBER

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
class AMBER(Engine):

    """ Engine for carrying out general purpose AMBER calculations. """

    def __init__(self, name="amber", **kwargs):
        ## Keyword args that aren't in this list are filtered out.
        self.valkwd = ['amberhome', 'pdb', 'mol2', 'frcmod', 'leapcmd']
        super(AMBER,self).__init__(name=name, **kwargs)

    def setopts(self, **kwargs):
        
        """ Called by __init__ ; Set AMBER-specific options. """

        ## The directory containing TINKER executables (e.g. dynamic)
        if 'amberhome' in kwargs:
            self.amberhome = kwargs['amberhome']
            if not os.path.exists(os.path.join(self.amberhome,"sander")):
                warn_press_key("The 'sander' executable indicated by %s doesn't exist! (Check amberhome)" \
                                   % os.path.join(self.amberhome,"sander"))
        else:
            warn_once("The 'amberhome' option was not specified; using default.")
            if which('sander') == '':
                warn_press_key("Please add AMBER executables to the PATH or specify amberhome.")
            self.amberhome = os.path.split(which('sander'))[0]
        
        with wopen('.quit.leap') as f:
            print >> f, 'quit'

        # AMBER search path
        self.spath = []
        for line in self.callamber('tleap -f .quit.leap'):
            if 'Adding' in line and 'to search path' in line:
                self.spath.append(line.split('Adding')[1].split()[0])
        os.remove('.quit.leap')

    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)

        # Write the PDB that AMBER is going to read in.
        # This may or may not be identical to the one used to initialize the engine.
        # self.mol.write('%s.pdb' % self.name)
        # self.abspdb = os.path.abspath('%s.pdb' % self.name)

    def callamber(self, command, stdin=None, print_to_screen=False, print_command=False, **kwargs):

        """ Call TINKER; prepend the amberhome to calling the TINKER program. """

        csplit = command.split()
        # Sometimes the engine changes dirs and the inpcrd/prmtop go missing, so we link it.
        # Prepend the AMBER path to the program call.
        prog = os.path.join(self.amberhome, "bin", csplit[0])
        csplit[0] = prog
        # No need to catch exceptions since failed AMBER calculations will return nonzero exit status.
        o = _exec(' '.join(csplit), stdin=stdin, print_to_screen=print_to_screen, print_command=print_command, rbytes=1024, **kwargs)
        return o

    def leap(self, name=None, delcheck=False):
        if not os.path.exists(self.leapcmd):
            LinkFile(self.absleap, self.leapcmd)
        pdb = os.path.basename(self.abspdb)
        if not os.path.exists(pdb):
            LinkFile(self.abspdb, pdb)
        if name == None: name = self.name
        write_leap(self.leapcmd, mol2=self.mol2, frcmod=self.frcmod, pdb=pdb, prefix=name, spath=self.spath, delcheck=delcheck)
        self.callamber("tleap -f %s_" % self.leapcmd)

    def prepare(self, pbc=False, **kwargs):

#.........这里部分代码省略.........
开发者ID:albaugh,项目名称:forcebalance,代码行数:103,代码来源:amberio.py

示例10: TINKER

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]

#.........这里部分代码省略.........
                warn_once("Deleting PBC options from the .key file.")
                tk_opts['a-axis'] = None
                tk_opts['b-axis'] = None
                tk_opts['c-axis'] = None
                tk_opts['alpha'] = None
                tk_opts['beta'] = None
                tk_opts['gamma'] = None
        if pbc:
            if (not keypbc) and 'boxes' not in self.mol.Data:
                logger.error("Periodic boundary conditions require either (1) a-axis to be in the .key file or (b) boxes to be in the coordinate file.\n")
                raise RuntimeError
        self.pbc = pbc
        if pbc:
            tk_opts['ewald'] = ''
            if minbox <= 10:
                warn_press_key("Periodic box is set to less than 10 Angstroms across")
            # TINKER likes to use up to 7.0 Angstrom for PME cutoffs
            rpme = 0.05*(float(int(minbox - 1))) if minbox <= 15 else 7.0
            tk_defs['ewald-cutoff'] = "%f" % rpme
            # TINKER likes to use up to 9.0 Angstrom for vdW cutoffs
            rvdw = 0.05*(float(int(minbox - 1))) if minbox <= 19 else 9.0
            tk_defs['vdw-cutoff'] = "%f" % rvdw
            if (minbox*0.5 - rpme) > 2.5 and (minbox*0.5 - rvdw) > 2.5:
                tk_defs['neighbor-list'] = ''
            elif (minbox*0.5 - rpme) > 2.5:
                tk_defs['mpole-list'] = ''
        else:
            tk_opts['ewald'] = None
            tk_opts['ewald-cutoff'] = None
            tk_opts['vdw-cutoff'] = None
            # This seems to have no effect on the kinetic energy.
            # tk_opts['remove-inertia'] = '0'

        write_key("%s.key" % self.name, tk_opts, os.path.join(self.srcdir, self.key) if self.key else None, tk_defs, verbose=False, prmfnm=prmfnm)
        self.abskey = os.path.abspath("%s.key")

        self.mol[0].write(os.path.join("%s.xyz" % self.name), ftype="tinker")

        ## If the coordinates do not come with TINKER suffixes then throw an error.
        self.mol.require('tinkersuf')

        ## Call analyze to read information needed to build the atom lists.
        o = self.calltinker("analyze %s.xyz P,C" % (self.name), stdin="ALL")

        ## Parse the output of analyze.
        mode = 0
        self.AtomMask = []
        self.AtomLists = defaultdict(list)
        ptype_dict = {'atom': 'A', 'vsite': 'D'}
        G = nx.Graph()
        for line in o:
            s = line.split()
            if len(s) == 0: continue
            if "Atom Type Definition Parameters" in line:
                mode = 1
            if mode == 1:
                if isint(s[0]): mode = 2
            if mode == 2:
                if isint(s[0]):
                    mass = float(s[5])
                    self.AtomLists['Mass'].append(mass)
                    if mass < 1.0:
                        # Particles with mass less than one count as virtual sites.
                        self.AtomLists['ParticleType'].append('D')
                    else:
                        self.AtomLists['ParticleType'].append('A')
开发者ID:albaugh,项目名称:forcebalance,代码行数:70,代码来源:tinkerio.py

示例11: len

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
                espval.append(float(s[3]))
            elif len(espxyz) > 0:
                # After reading in a block of ESPs, don't read any more.
                ESPMode = -1 
        if line.strip().startswith("Geometry (in Angstrom)"):
            XMode = 1
            EMode = len(elem) == 0
        if 'Electrostatic Potential' in line.strip() and ESPMode == 0:
            ESPMode = 1
    if len(xyzs) == 0:
        raise Exception('%s has length zero' % psiout)
    return xyzs, elem, espxyz, espval

xyzs, elem, espxyz, espval = read_psi_xyzesp(sys.argv[1])

M = Molecule()
M.xyzs = xyzs
M.elem = elem
M.write('%s.xyz' % os.path.splitext(sys.argv[1])[0])

EM = Molecule()
EM.xyzs = [np.array(espxyz) * 0.52917721092]
EM.elem = ['H' for i in range(len(espxyz))]
EM.write('%s.espx' % os.path.splitext(sys.argv[1])[0], ftype="xyz")

M.qm_espxyzs = EM.xyzs
M.qm_espvals = [np.array(espval)]
M.write("qdata.txt")

np.savetxt('%s.esp' % os.path.splitext(sys.argv[1])[0], espval)
开发者ID:leeping,项目名称:forcebalance,代码行数:32,代码来源:read-esp.py

示例12: main

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
def main():
    topfnm = argv[3] if len(argv) >= 4 else None
    M = Molecule(argv[1], top=topfnm)
    M.write(argv[2])
开发者ID:falagara,项目名称:forcebalance,代码行数:6,代码来源:filecnv.py

示例13: int

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
# Psi4 output file.
psiout = sys.argv[1]

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

frqs, modes, elem, xyz = read_frq_psi(psiout)

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)

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))
]

M.write(os.path.splitext(psiout)[0] + ".mode%03i.xyz" % modenum)
开发者ID:mesoniat,项目名称:forcebalance,代码行数:31,代码来源:anifrq.py

示例14: int

# 需要导入模块: from forcebalance.molecule import Molecule [as 别名]
# 或者: from forcebalance.molecule.Molecule import write [as 别名]
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:
                Cnvgd[s[0]] = float(s[-3])
开发者ID:ChayaSt,项目名称:open-forcefield-group,代码行数:32,代码来源:opt-qchem.py


注:本文中的forcebalance.molecule.Molecule.write方法示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。