Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Extract the chemical element.
Pelem.append(pid.split(':')[1].split(',')[0].split('=')[1])
Pelem = set(Pelem)
if len(self.Elements.intersection(Pelem)) == 0:
self.call_derivatives[p] = False
## Psi4 basis set file
gbslist = [i for i in self.FF.fnms if os.path.splitext(i)[1] == '.gbs']
if len(gbslist) != 1:
warn_press_key("In %s, you should only have exactly one .gbs file in the list of force field files!" % __file__)
self.GBSfnm = gbslist[0]
## Psi4 input file for calculation of linear dependencies
## This is actually a file in 'forcefield' until we can figure out a better system
if CheckBasis():
datlist = [i for i in self.FF.fnms if os.path.splitext(i)[1] == '.dat']
if len(datlist) != 1:
warn_press_key("In %s, you should only have exactly one .dat file in the list of force field files!" % __file__)
self.DATfnm = datlist[0]
## Prepare the temporary directory
self.prepare_temp_directory(options,tgt_opts)
def make_redirect(self,mvals):
Groups = defaultdict(list)
for p, pid in enumerate(self.plist):
if 'Exponent' not in pid or len(pid.split()) != 1:
warn_press_key("Fusion penalty currently implemented only for basis set optimizations, where parameters are like this: Exponent:Elem=H,AMom=D,Bas=0,Con=0")
Data = dict([(i.split('=')[0],i.split('=')[1]) for i in pid.split(':')[1].split(',')])
if 'Con' not in Data or Data['Con'] != '0':
warn_press_key("More than one contraction coefficient found! You should expect the unexpected")
key = Data['Elem']+'_'+Data['AMom']
Groups[key].append(p)
#pvals = self.FF.create_pvals(mvals)
#print "pvals: ", pvals
pvals = self.create_pvals(mvals)
logger.info("pvals:\n")
logger.info(str(pvals) + '\n')
Thresh = 1e-4
for gnm, pidx in Groups.items():
# The group of parameters for a particular element / angular momentum.
pvals_grp = pvals[pidx]
# The order that the parameters come in.
Order = argsort(pvals_grp)
def setopts(self, **kwargs):
""" Called by __init__ ; Set AMBER-specific options. """
if 'amberhome' in kwargs:
self.amberhome = kwargs['amberhome']
if not os.path.exists(os.path.join(self.amberhome, "bin", "sander")):
warn_press_key("The 'sander' executable indicated by %s doesn't exist! (Check amberhome)" \
% os.path.join(self.amberhome,"bin","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]
self.have_pmemd_cuda = False
if os.path.exists(os.path.join(self.amberhome, "bin", "pmemd.cuda")):
self.callamber('pmemd.cuda -h', persist=True)
if _exec.returncode != 0:
warn_press_key("pmemd.cuda gave a nonzero returncode; CUDA environment variables likely missing")
else:
logger.info("pmemd.cuda is available, using CUDA to accelerate calculations.\n")
self.have_pmemd_cuda = True
parameters like C_ , G_ .. or simply unit conversions, you get the idea.
@warning I don't think the multiplier actually works for analytic derivatives unless
the interaction calculator knows the multiplier as well. I'm sure I can make this
work in the future if necessary.
@param[in] ffname Name of the force field file
"""
fftype = determine_fftype(ffname)
ffname = ffname.split(':')[0]
# Set the Tinker PRM file, which will help for running programs like "analyze".
if fftype == "tinker":
if hasattr(self, "tinkerprm"):
warn_press_key("There should only be one TINKER parameter file")
else:
self.tinkerprm = ffname
# Set the OpenMM XML file, which will help for running OpenMM.
if fftype == "openmm":
if hasattr(self, "openmmxml"):
warn_press_key("There should only be one OpenMM XML file - confused!!")
else:
self.openmmxml = ffname
# Determine the appropriate parser from the FF_IOModules dictionary.
# If we can't figure it out, then use the base reader, it ain't so bad. :)
Reader = FF_IOModules.get(fftype, forcebalance.BaseReader)
# Open the force field using an absolute path and read its contents into memory.
absff = os.path.join(self.root,self.ffdir,ffname)
def find_spacings(self):
Groups = defaultdict(list)
for p, pid in enumerate(self.plist):
if 'Exponent' not in pid or len(pid.split()) != 1:
warn_press_key("Fusion penalty currently implemented only for basis set optimizations, where parameters are like this: Exponent:Elem=H,AMom=D,Bas=0,Con=0")
Data = dict([(i.split('=')[0],i.split('=')[1]) for i in pid.split(':')[1].split(',')])
if 'Con' not in Data or Data['Con'] != '0':
warn_press_key("More than one contraction coefficient found! You should expect the unexpected")
key = Data['Elem']+'_'+Data['AMom']
Groups[key].append(p)
pvals = self.create_pvals(zeros(self.np,dtype=float))
logger.info("pvals:\n")
logger.info(str(pvals) + '\n')
spacdict = {}
for gnm, pidx in Groups.items():
# The group of parameters for a particular element / angular momentum.
pvals_grp = pvals[pidx]
# The order that the parameters come in.
Order = argsort(pvals_grp)
spacs = []
for p in range(len(Order) - 1):
# The pointers to the parameter indices.
Emin = min(Eig)
if Emin < self.eps: # Mix in SD step if Hessian minimum eigenvalue is negative
# Experiment.
Adj = max(self.eps, 0.01*abs(Emin)) - Emin
logger.info("Hessian has a small or negative eigenvalue (%.1e), mixing in some steepest descent (%.1e) to correct this.\n" % (Emin, Adj))
logger.info("Eigenvalues are:\n")
pvec1d(Eig)
H += Adj*np.eye(H.shape[0])
if self.bhyp:
G = np.delete(G, self.excision)
H = np.delete(H, self.excision, axis=0)
H = np.delete(H, self.excision, axis=1)
xkd = np.delete(xk, self.excision)
if self.Objective.Penalty.fmul != 0.0:
warn_press_key("Using the multiplicative hyperbolic penalty is discouraged!")
# This is the gradient and Hessian without the contributions from the hyperbolic constraint.
Obj0 = {'X':X,'G':G,'H':H}
class Hyper(object):
def __init__(self, HL, Penalty):
self.H = HL.copy()
self.dx = 1e10 * np.ones(len(HL))
self.Val = 0
self.Grad = np.zeros(len(HL))
self.Hess = np.zeros((len(HL),len(HL)))
self.Penalty = Penalty
self.iter = 0
def _compute(self, dx):
self.dx = dx.copy()
Tmp = np.dot(self.H, dx)
Reg_Term = self.Penalty.compute(xkd+flat(dx), Obj0)
self.Val = (X + np.dot(dx, G) + 0.5*np.dot(dx,Tmp) + Reg_Term[0] - data['X'])
def prepare_temp_directory(self, options, tgt_opts):
os.environ["GMX_NO_SOLV_OPT"] = "TRUE"
os.environ["GMX_NO_ALLVSALL"] = "TRUE"
abstempdir = os.path.join(self.root,self.tempdir)
if options['gmxpath'] == None or options['gmxsuffix'] == None:
warn_press_key('Please set the options gmxpath and gmxsuffix in the input file!')
if not os.path.exists(os.path.join(options['gmxpath'],"mdrun"+options['gmxsuffix'])):
warn_press_key('The mdrun executable pointed to by %s doesn\'t exist! (Check gmxpath and gmxsuffix)' % os.path.join(options['gmxpath'],"mdrun"+options['gmxsuffix']))
# Link the necessary programs into the temporary directory
LinkFile(os.path.join(options['gmxpath'],"mdrun"+options['gmxsuffix']),os.path.join(abstempdir,"mdrun"))
LinkFile(os.path.join(options['gmxpath'],"grompp"+options['gmxsuffix']),os.path.join(abstempdir,"grompp"))
LinkFile(os.path.join(options['gmxpath'],"trjconv"+options['gmxsuffix']),os.path.join(abstempdir,"trjconv"))
# Link the run files
LinkFile(os.path.join(self.root,self.tgtdir,"conf.gro"),os.path.join(abstempdir,"conf.gro"))
LinkFile(os.path.join(self.root,self.tgtdir,"grompp.mdp"),os.path.join(abstempdir,"grompp.mdp"))
LinkFile(os.path.join(self.root,self.tgtdir,"topol.top"),os.path.join(abstempdir,"topol.top"))
Emin = min(Eig)
if Emin < self.eps: # Mix in SD step if Hessian minimum eigenvalue is negative
# Experiment.
Adj = max(self.eps, 0.01*abs(Emin)) - Emin
logger.info("Hessian has a small or negative eigenvalue (%.1e), mixing in some steepest descent (%.1e) to correct this.\n" % (Emin, Adj))
logger.info("Eigenvalues are:\n") ###
pvec1d(Eig) ###
H += Adj*np.eye(H.shape[0])
if self.bhyp:
G = np.delete(G, self.excision)
H = np.delete(H, self.excision, axis=0)
H = np.delete(H, self.excision, axis=1)
xkd = np.delete(xk, self.excision)
if self.Objective.Penalty.fmul != 0.0:
warn_press_key("Using the multiplicative hyperbolic penalty is discouraged!")
# This is the gradient and Hessian without the contributions from the hyperbolic constraint.
Obj0 = {'X':X,'G':G,'H':H}
class Hyper(object):
def __init__(self, HL, Penalty):
self.H = HL.copy()
self.dx = 1e10 * np.ones(len(HL),dtype=float)
self.Val = 0
self.Grad = np.zeros(len(HL),dtype=float)
self.Hess = np.zeros((len(HL),len(HL)),dtype=float)
self.Penalty = Penalty
def _compute(self, dx):
self.dx = dx.copy()
Tmp = np.mat(self.H)*col(dx)
Reg_Term = self.Penalty.compute(xkd+flat(dx), Obj0)
self.Val = (X + np.dot(dx, G) + 0.5*row(dx)*Tmp + Reg_Term[0] - data['X'])[0,0]
self.Grad = flat(col(G) + Tmp) + Reg_Term[1]
elif self.adapt_fac != 0.0:
detail = "(Adaptive Trust Radius)"
else:
detail = "(Trust Radius)"
printcool("Main Optimizer \n%s Method %s\n\n"
"\x1b[0mConvergence criteria (%i of 3 needed):\n"
"\x1b[0mObjective Function : %.3e\n"
"\x1b[0mNorm of Gradient : %.3e\n"
"\x1b[0mParameter step size : %.3e" %
("BFGS" if b_BFGS else "Newton-Raphson", detail, self.criteria,
self.convergence_objective, self.convergence_gradient,
self.convergence_step), ansi=1, bold=1)
# Print a warning if optimization is unlikely to converge
if self.uncert and self.convergence_objective < 1e-3:
warn_press_key("Condensed phase targets detected - may not converge with current choice of"
" convergence_objective (%.e)\nRecommended range is 1e-2 - 1e-1 for this option." % self.convergence_objective)
#========================#
#| Initialize variables |#
#========================#
# Order of derivatives
Ord = 1 if b_BFGS else 2
# Iteration number counter.
global ITERATION
ITERATION = self.iterinit
self.iteration = self.iterinit
# Indicates if the optimization step was "good" (i.e. not rejected).
self.set_goodstep(1)
# Indicates if the optimization is currently at the lowest value of the objective function so far.
Best_Step = 1
# Objective function history.
def build_pid(self, pfld):
if self.this_section not in ["global", "atom", "bond", "offdiag", "angle", "torsion", "hbond"]:
logger.info("%s\n" % re.sub(" *\n$", "", self.line))
nif.warn_press_key("This line is not part of a section that contains parameters.")
SectionID = self.this_section.capitalize()
AtomsID = '-'.join(self.involved)
if pfld >= len(self.field_names):
logger.info("%s\n" % re.sub(" *\n$", "", self.line))
nif.warn_press_key("Field number %i is not recognized in the line above, will lead to an error.\nThe fields for this line are: %s" % (pfld, str(self.field_names)))
ParamID = self.field_names[pfld]
if ParamID in UnusedParamNames[self.this_section]:
logger.warning("Field number %i has parameter type %s which is not supposed to be parameterized\n" % (pfld, ParamID))
print SectionID, AtomsID, ParamID
if self.this_section == "global":
return "/".join([SectionID, ParamID])
else:
return "/".join([SectionID, ParamID, AtomsID])