How to use the forcebalance.molecule.Molecule function in forcebalance

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:

def multipole_moments(self, shot=0, optimize=True, polarizability=False):
        
        """ Return the multipole moments of the 1st snapshot in Debye and Buckingham units. """
        
        if not self.double:
            warn_once("Single-precision GROMACS detected - recommend that you use double precision build.")

        if polarizability:
            raise NotImplementedError

        if optimize: 
            self.optimize(shot)
            M = Molecule("%s-min.gro" % self.name)
        else:
            self.mol[shot].write("%s.gro" % self.name)
            M = Molecule("%s.gro" % self.name)
        
        #-----
        # g_dipoles uses a different reference point compared to TINKER
        #-----
        # self.callgmx("g_dipoles -s %s-d.tpr -f %s-d.gro -o %s-d.xvg -xvg no" % (self.name, self.name, self.name), stdin="System\n")
        # Dips = np.array([[float(i) for i in line.split()[1:4]] for line in open("%s-d.xvg" % self.name)])[0]
        #-----
        
        ea_debye = 4.803204255928332 # Conversion factor from e*nm to Debye
        q = self.get_charges()
        x = M.xyzs[0] - M.center_of_mass()[0]

        xx, xy, xz, yy, yz, zz = (x[:,i]*x[:,j] for i, j in [(0,0),(0,1),(0,2),(1,1),(1,2),(2,2)])

if self.have_constraints:
                algorithm = "steep"
            else:
                algorithm = "l-bfgs"
            # Arguments for running minimization.
            min_opts = {"integrator" : algorithm, "emtol" : crit, "nstxout" : 0, "nstfout" : 0, "nsteps" : 10000, "nstenergy" : 1}

        edit_mdp(fin="%s.mdp" % self.name, fout="%s-min.mdp" % self.name, options=min_opts)

        self.warngmx("grompp -c %s.gro -p %s.top -f %s-min.mdp -o %s-min.tpr" % (self.name, self.name, self.name, self.name))
        self.callgmx("mdrun -deffnm %s-min -nt 1" % self.name)
        self.callgmx("trjconv -f %s-min.trr -s %s-min.tpr -o %s-min.gro -ndec 9" % (self.name, self.name, self.name), stdin="System")
        self.callgmx("g_energy -xvg no -f %s-min.edr -o %s-min-e.xvg" % (self.name, self.name), stdin='Potential')
        
        E = float(open("%s-min-e.xvg" % self.name).readlines()[-1].split()[1])
        M = Molecule("%s.gro" % self.name, build_topology=False) + Molecule("%s-min.gro" % self.name, build_topology=False)
        if not self.pbc:
            M.align(center=False)
        rmsd = M.ref_rmsd(0)[1]
        M[1].write("%s-min.gro" % self.name)

        return E / 4.184, rmsd, M[1]

# Figure out where the trajectory information is coming from
        # First priority: Passed as input to trajfnm
        # Second priority: Using self.trajectory filename attribute
        # Third priority: Using internal Molecule object
        # 0 = using internal Molecule object
        # 1 = using NetCDF trajectory format
        mode = 0
        if traj_fnm is None and hasattr(self, 'trajectory'):
            traj_fnm = self.trajectory
        if traj_fnm is not None:
            try:
                nc = netcdf_file(traj_fnm, 'r')
                mode = 1
            except TypeError:
                print("Failed to load traj as netcdf, trying to load as Molecule object")
                mol = Molecule(traj_fnm)
        else:
            mol = self.mol

        def get_coord_box(i):
            box = None
            if mode == 0:
                coords = mol.xyzs[i]
                if self.pbc:
                    box = [mol.boxes[i].a, mol.boxes[i].b, mol.boxes[i].c,
                           mol.boxes[i].alpha, mol.boxes[i].beta, mol.boxes[i].gamma]
            elif mode == 1:
                coords = nc.variables['coordinates'].data[i].copy()
                if self.pbc:
                    cl = nc.variables['cell_lengths'].data[i].copy()
                    ca = nc.variables['cell_angles'].data[i].copy()
                    box = [cl[0], cl[1], cl[2], ca[0], ca[1], ca[2]]

#!/usr/bin/env python

from __future__ import division
from __future__ import print_function
from builtins import range
import os, sys, re
import numpy as np
from forcebalance.molecule import Molecule
from forcebalance.readfrq import read_frq_tc

# TeraChem frequency output file.
tcout = sys.argv[1]

# Starting coordinate file.
M = Molecule(sys.argv[2])
xyz = M.xyzs[0]
M.xyzs = []
M.comms = []

# Mode number, starting from 1.
modenum = int(sys.argv[3])
if modenum == 0:
    raise RuntimeError('Mode numbers start from 1')

frqs, modes = read_frq_tc(tcout)

xmode = modes[modenum - 1].reshape(-1,3)

xmodeNorm = np.array([np.linalg.norm(i) for i in xmode])
idxMax = np.argmax(xmodeNorm)
print("In mode #%i, the largest displacement comes from atom #%i (%s); norm %.3f" % (modenum, idxMax+1, M.elem[idxMax], np.max(xmodeNorm)))

def readsrc(self, **kwargs):
        """ Called by __init__ ; read files from the source directory. """

        ## Attempt to determine file names of .gro, .top, and .mdp files
        self.top = onefile(kwargs.get('gmx_top'), 'top')
        self.mdp = onefile(kwargs.get('gmx_mdp'), 'mdp')
        if 'mol' in kwargs:
            self.mol = kwargs['mol']
        else:
            self.mol = Molecule(onefile(kwargs.get('coords'), 'gro', err=True))

# 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)
                # LinkFile(os.path.join(self.root,self.tgtdir,sysopt['keyfile']),os.path.join(abstempdir,sysname+".key"))

def readsrc(self, **kwargs):

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

        ## Attempt to determine file names of .gro, .top, and .mdp files
        self.top = onefile('top', kwargs['gmx_top'] if 'gmx_top' in kwargs else None)
        self.mdp = onefile('mdp', kwargs['gmx_mdp'] if 'gmx_mdp' in kwargs else None)
        if 'mol' in kwargs:
            self.mol = kwargs['mol']
        elif 'coords' in kwargs and os.path.exists(kwargs['coords']):
            self.mol = Molecule(kwargs['coords'])
        else:
            grofile = onefile('gro')
            self.mol = Molecule(grofile)

minimize = TgtOptions.get('minimize_energy', 1)

    # Print all options.
    printcool_dictionary(TgtOptions, title="Options from ForceBalance")
    lipid_snapshots = int((lipid_nsteps * lipid_timestep / 1000) / lipid_intvl)
    lipid_iframes = int(1000 * lipid_intvl / lipid_timestep)
    logger.info("For the condensed phase system, I will collect %i snapshots spaced apart by %i x %.3f fs time steps\n" \
        % (lipid_snapshots, lipid_iframes, lipid_timestep))
    if lipid_snapshots < 2:
        raise Exception('Please set the number of lipid time steps so that you collect at least two snapshots (minimum %i)' \
                            % (2000 * int(lipid_intvl/lipid_timestep)))

    #----
    # Loading coordinates
    #----
    ML = Molecule(lipid_fnm, toppbc=True)
    # Determine the number of molecules in the condensed phase coordinate file.
    NMol = len(ML.molecules)

    #----
    # Setting up MD simulations
    #----
    EngOpts = OrderedDict()
    EngOpts["lipid"] = OrderedDict([("coords", lipid_fnm), ("mol", ML), ("pbc", True)])
    if "nonbonded_cutoff" in TgtOptions:
        EngOpts["lipid"]["nonbonded_cutoff"] = TgtOptions["nonbonded_cutoff"]
    if "vdw_cutoff" in TgtOptions:
        EngOpts["lipid"]["vdw_cutoff"] = TgtOptions["vdw_cutoff"]
    GenOpts = OrderedDict([('FF', FF)])
    if engname == "openmm":
        # OpenMM-specific options
        EngOpts["liquid"]["platname"] = TgtOptions.get("platname", 'CUDA')

def read_quantum():
        Result = None
        os.chdir('calcs')
        for i in range(M.ns):
            dnm = eval(formstr % i)
            print "\rNow in directory %i" % i,
            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