Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
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