Api
pygsm.growing_string_methods.de_gsm
🔗
Classes🔗
DE_GSM
🔗
Attributes🔗
dEs
inherited
property
readonly
🔗
dqmaga
inherited
property
writable
🔗
emax
inherited
property
readonly
🔗
energies
inherited
property
readonly
🔗
Energies of string
geometries
inherited
property
readonly
🔗
gradrmss
inherited
property
readonly
🔗
ictan
inherited
property
writable
🔗
npeaks
inherited
property
readonly
🔗
TSnode
inherited
property
readonly
🔗
The current node with maximum energy
Methods🔗
__init__(self, options)
special
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def __init__(
self,
options,
):
super(DE_GSM,self).__init__(options)
print(" Assuming primitives are union!")
print(" number of primitives is", self.nodes[0].num_primitives)
add_GSM_nodeP(self, newnodes = 1)
inherited
🔗
Add a node between endpoints on the product side, should only be called inside GSM
Source code in pygsm/growing_string_methods/de_gsm.py
def add_GSM_nodeP(self,newnodes=1):
'''
Add a node between endpoints on the product side, should only be called inside GSM
'''
printcool("Adding product node")
if self.current_nnodes+newnodes > self.nnodes:
raise ValueError("Adding too many nodes, cannot interpolate")
for i in range(newnodes):
#self.nodes[-self.nP-1] = BaseClass.add_node(self.nnodes-self.nP,self.nnodes-self.nP-1,self.nnodes-self.nP)
n1=self.nnodes-self.nP
n2=self.nnodes-self.nP-1
n3=self.nR-1
print(" adding node: %i between %i %i from %i" %(n2,n1,n3,n1))
if self.nnodes - self.current_nnodes > 1:
stepsize = 1./float(self.nnodes-self.current_nnodes+1)
else:
stepsize = 0.5
self.nodes[-self.nP-1] = GSM.add_node(
self.nodes[n1],
self.nodes[n3],
stepsize,
n2
)
if self.nodes[-self.nP-1]==None:
raise Exception('Ran out of space')
self.optimizer[n2].DMAX = self.optimizer[n1].DMAX
self.current_nnodes+=1
self.nP+=1
print(" nn=%i,nP=%i" %(self.current_nnodes,self.nP))
self.active[-self.nP] = True
# align center of mass and rotation
#print("%i %i %i" %(n1,n3,n2))
#print(" Aligning")
self.nodes[-self.nP].xyz = self.com_rotate_move(n1,n3,n2)
#print(" getting energy for node %d: %5.4f" %(self.nnodes-self.nP,self.nodes[-self.nP].energy - self.nodes[0].V0))
return
add_GSM_nodeR(self, newnodes = 1)
inherited
🔗
Add a node between endpoints on the reactant side, should only be called inside GSM
Source code in pygsm/growing_string_methods/de_gsm.py
def add_GSM_nodeR(self,newnodes=1):
'''
Add a node between endpoints on the reactant side, should only be called inside GSM
'''
printcool("Adding reactant node")
if self.current_nnodes+newnodes > self.nnodes:
raise ValueError("Adding too many nodes, cannot interpolate")
for i in range(newnodes):
iR = self.nR-1
iP = self.nnodes-self.nP
iN = self.nR
print(" adding node: %i between %i %i from %i" %(iN,iR,iP,iR))
if self.nnodes - self.current_nnodes > 1:
stepsize = 1./float(self.nnodes-self.current_nnodes+1)
else:
stepsize = 0.5
self.nodes[self.nR] = GSM.add_node(
self.nodes[iR],
self.nodes[iP],
stepsize,
iN,
DQMAG_MAX = self.DQMAG_MAX,
DQMAG_MIN = self.DQMAG_MIN,
driving_coords = self.driving_coords,
)
if self.nodes[self.nR]==None:
raise Exception('Ran out of space')
if self.__class__.__name__!="DE_GSM":
ictan,bdist = self.get_tangent(
self.nodes[self.nR],
None,
driving_coords=self.driving_coords,
)
self.nodes[self.nR].bdist = bdist
self.optimizer[self.nR].DMAX = self.optimizer[self.nR-1].DMAX
self.current_nnodes+=1
self.nR+=1
print(" nn=%i,nR=%i" %(self.current_nnodes,self.nR))
self.active[self.nR-1] = True
# align center of mass and rotation
#print("%i %i %i" %(iR,iP,iN))
#print(" Aligning")
self.nodes[self.nR-1].xyz = self.com_rotate_move(iR,iP,iN)
add_GSM_nodes(self, newnodes = 1)
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def add_GSM_nodes(self,newnodes=1):
if self.current_nnodes+newnodes > self.nnodes:
print("Adding too many nodes, cannot add_GSM_node")
sign = -1
for i in range(newnodes):
sign *= -1
if sign == 1:
self.add_GSM_nodeR()
else:
self.add_GSM_nodeP()
add_node(nodeR, nodeP, stepsize, node_id, **kwargs)
inherited
🔗
Add a node between nodeR and nodeP or if nodeP is none use driving coordinate to add new node
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def add_node(
nodeR,
nodeP,
stepsize,
node_id,
**kwargs
):
'''
Add a node between nodeR and nodeP or if nodeP is none use driving coordinate to add new node
'''
#get driving coord
driving_coords = kwargs.get('driving_coords',None)
DQMAG_MAX =kwargs.get('DQMAG_MAX',0.8)
DQMAG_MIN =kwargs.get('DQMAG_MIN',0.2)
if nodeP is None:
if driving_coords is None:
raise RuntimeError("You didn't supply a driving coordinate and product node is None!")
BDISTMIN=0.05
ictan,bdist = GSM.get_tangent(nodeR,None,driving_coords=driving_coords)
if bdist<BDISTMIN:
print("bdist too small %.3f" % bdist)
return None
new_node = Molecule.copy_from_options(nodeR,new_node_id=node_id)
Vecs = new_node.update_coordinate_basis(constraints=ictan)
constraint = new_node.constraints[:,0]
sign=-1.
dqmag_scale=1.5
minmax = DQMAG_MAX - DQMAG_MIN
a = bdist/dqmag_scale
if a>1.:
a=1.
dqmag = sign*(DQMAG_MIN+minmax*a)
if dqmag > DQMAG_MAX:
dqmag = DQMAG_MAX
print(" dqmag: %4.3f from bdist: %4.3f" %(dqmag,bdist))
dq0 = dqmag*constraint
print(" dq0[constraint]: %1.3f" % dqmag)
new_node.update_xyz(dq0)
new_node.bdist = bdist
else:
ictan,_ = GSM.get_tangent(nodeR,nodeP)
Vecs = nodeR.update_coordinate_basis(constraints=ictan)
constraint = nodeR.constraints[:,0]
dqmag = np.linalg.norm(ictan)
print(" dqmag: %1.3f"%dqmag)
#sign=-1
sign=1.
dqmag *= (sign*stepsize)
print(" scaled dqmag: %1.3f"%dqmag)
dq0 = dqmag*constraint
old_xyz = nodeR.xyz.copy()
new_xyz = nodeR.coord_obj.newCartesian(old_xyz,dq0)
new_node = Molecule.copy_from_options(MoleculeA=nodeR,xyz=new_xyz,new_node_id=node_id)
return new_node
calc_optimization_metrics(nodes)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def calc_optimization_metrics(nodes):
'''
'''
nnodes = len(nodes)
rn3m6 = np.sqrt(3*nodes[0].natoms-6)
totalgrad = 0.0
gradrms = 0.0
sum_gradrms = 0.0
for i,ico in enumerate(nodes[1:nnodes-1]):
if ico!=None:
print(" node: {:02d} gradrms: {:.6f}".format(i,float(ico.gradrms)),end='')
if i%5 == 0:
print()
totalgrad += ico.gradrms*rn3m6
gradrms += ico.gradrms*ico.gradrms
sum_gradrms += ico.gradrms
print('')
#TODO wrong for growth
gradrms = np.sqrt(gradrms/(nnodes-2))
return totalgrad,gradrms,sum_gradrms
check_if_grown(self)
🔗
Check if the string is grown Returns True if grown
Source code in pygsm/growing_string_methods/de_gsm.py
def check_if_grown(self):
'''
Check if the string is grown
Returns True if grown
'''
return self.current_nnodes==self.nnodes
check_tscontinue(self)
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def check_tscontinue(self):
tol1=0.5
tol2=2.
allup=True
alldown=True
diss=False
energies = self.energies
nnodes = self.nnodes
# if allup
for n in range(1,nnodes):
if energies[n]+tol1<energies[n-1]:
allup=False
break
# alldown
for n in range(1,nnodes-1):
if energies[n+1]+tol1>energies[n]:
alldown=False
break
if allup or alldown:
self.tscontinue=False
return
com_rotate_move(self, iR, iP, iN)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def com_rotate_move(self,iR,iP,iN):
print(" aligning com and to Eckart Condition")
mfrac = 0.5
if self.nnodes - self.current_nnodes+1 != 1:
mfrac = 1./(self.nnodes - self.current_nnodes+1)
#if self.__class__.__name__ != "DE_GSM":
# # no "product" structure exists, use initial structure
# iP = 0
xyz0 = self.nodes[iR].xyz.copy()
xyz1 = self.nodes[iN].xyz.copy()
com0 = self.nodes[iR].center_of_mass
com1 = self.nodes[iN].center_of_mass
masses = self.nodes[iR].mass_amu
# From the old GSM code doesn't work
#com1 = mfrac*(com2-com0)
#print("com1")
#print(com1)
## align centers of mass
#xyz1 += com1
#Eckart_align(xyz1,xyz2,masses,mfrac)
# rotate to be in maximal coincidence with 0
# assumes iP i.e. 2 is also in maximal coincidence
U = rotate.get_rot(xyz0,xyz1)
xyz1 = np.dot(xyz1,U)
## align
#if self.nodes[iP] != None:
# xyz2 = self.nodes[iP].xyz.copy()
# com2 = self.nodes[iP].center_of_mass
# if abs(iN-iR) > abs(iN-iP):
# avg_com = mfrac*com2 + (1.-mfrac)*com0
# else:
# avg_com = mfrac*com0 + (1.-mfrac)*com2
# dist = avg_com - com1 #final minus initial
#else:
# dist = com0 - com1 #final minus initial
#print("aligning to com")
#print(dist)
#xyz1 += dist
return xyz1
default_options()
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def default_options():
if hasattr(GSM, '_default_options'): return GSM._default_options.copy()
opt = options.Options()
opt.add_option(
key='reactant',
required=True,
#allowed_types=[Molecule,wrappers.Molecule],
doc='Molecule object as the initial reactant structure')
opt.add_option(
key='product',
required=False,
#allowed_types=[Molecule,wrappers.Molecule],
doc='Molecule object for the product structure (not required for single-ended methods.')
opt.add_option(
key='nnodes',
required=False,
value=1,
allowed_types=[int],
doc="number of string nodes"
)
opt.add_option(
key='optimizer',
required=True,
doc='Optimzer object to use e.g. eigenvector_follow, conjugate_gradient,etc. \
most of the default options are okay for here since GSM will change them anyway',
)
opt.add_option(
key='driving_coords',
required=False,
value=[],
allowed_types=[list],
doc='Provide a list of tuples to select coordinates to modify atoms\
indexed at 1')
opt.add_option(
key='CONV_TOL',
value=0.0005,
required=False,
allowed_types=[float],
doc='Convergence threshold'
)
opt.add_option(
key='CONV_gmax',
value=0.0005,
required=False,
allowed_types=[float],
doc='Convergence threshold'
)
opt.add_option(
key='CONV_Ediff',
value=0.1,
required=False,
allowed_types=[float],
doc='Convergence threshold'
)
opt.add_option(
key='CONV_dE',
value=0.5,
required=False,
allowed_types=[float],
doc='Convergence threshold'
)
opt.add_option(
key='ADD_NODE_TOL',
value=0.1,
required=False,
allowed_types=[float],
doc='Convergence threshold')
opt.add_option(
key="growth_direction",
value=0,
required=False,
doc="how to grow string,0=Normal,1=from reactant"
)
opt.add_option(
key="DQMAG_MAX",
value=0.8,
required=False,
doc="max step along tangent direction for SSM"
)
opt.add_option(
key="DQMAG_MIN",
value=0.2,
required=False,
doc=""
)
opt.add_option(
key='print_level',
value=1,
required=False
)
opt.add_option(
key='use_multiprocessing',
value=False,
doc='Use python multiprocessing module, an OpenMP like implementation \
that parallelizes optimization cycles on a single compute node'
)
opt.add_option(
key="BDIST_RATIO",
value=0.5,
required=False,
doc="SE-Crossing uses this \
bdist must be less than 1-BDIST_RATIO of initial bdist in order to be \
to be considered grown.",
)
opt.add_option(
key='ID',
value=0,
required=False,
doc='A number for identification of Strings'
)
opt.add_option(
key='interp_method',
value='DLC',
allowed_values=['Geodesic','DLC'],
required=False,
doc='Which reparameterization method to use',
)
opt.add_option(
key='noise',
value=1.0,
allowed_types=[float],
required=False,
doc='Noise to check for intermediate',
)
GSM._default_options = opt
return GSM._default_options.copy()
find_peaks(self, rtype)
inherited
🔗
This doesnt actually calculate peaks, it calculates some other thing
Source code in pygsm/growing_string_methods/de_gsm.py
def find_peaks(self,rtype):
'''
This doesnt actually calculate peaks, it calculates some other thing
'''
#rtype 1: growing
#rtype 2: opting
#rtype 3: intermediate check
if rtype==1:
nnodes=self.nR
elif rtype==2 or rtype==3:
nnodes=self.nnodes
else:
raise ValueError("find peaks bad input")
#if rtype==1 or rtype==2:
# print "Energy"
alluptol=0.1
alluptol2=0.5
allup=True
diss=False
energies = self.energies
for n in range(1,len(energies[:nnodes])):
if energies[n]+alluptol<energies[n-1]:
allup=False
break
if energies[nnodes-1]>15.0:
if nnodes-3>0:
if ((energies[nnodes-1]-energies[nnodes-2])<alluptol2 and
(energies[nnodes-2]-energies[nnodes-3])<alluptol2 and
(energies[nnodes-3]-energies[nnodes-4])<alluptol2):
print(" possible dissociative profile")
diss=True
print(" nnodes ",nnodes)
print(" all uphill? ",allup)
print(" dissociative? ",diss)
npeaks1=0
npeaks2=0
minnodes=[]
maxnodes=[]
if energies[1]>energies[0]:
minnodes.append(0)
if energies[nnodes-1]<energies[nnodes-2]:
minnodes.append(nnodes-1)
for n in range(self.n0,nnodes-1):
if energies[n+1]>energies[n]:
if energies[n]<energies[n-1]:
minnodes.append(n)
if energies[n+1]<energies[n]:
if energies[n]>energies[n-1]:
maxnodes.append(n)
print(" min nodes ",minnodes)
print(" max nodes ", maxnodes)
npeaks1 = len(maxnodes)
#print "number of peaks is ",npeaks1
ediff=0.5
PEAK4_EDIFF = 2.0
if rtype==1:
ediff=1.
if rtype==3:
ediff=PEAK4_EDIFF
if rtype==1:
nmax = np.argmax(energies[:self.nR])
emax = float(max(energies[:self.nR]))
else:
emax = float(max(energies))
nmax = np.argmax(energies)
print(" emax and nmax in find peaks %3.4f,%i " % (emax,nmax))
#check if any node after peak is less than 2 kcal below
for n in maxnodes:
diffs=( energies[n]-e>ediff for e in energies[n:nnodes])
if any(diffs):
found=n
npeaks2+=1
npeaks = npeaks2
print(" found %i significant peak(s) TOL %3.2f" %(npeaks,ediff))
#handle dissociative case
if rtype==3 and npeaks==1:
nextmin=0
for n in range(found,nnodes-1):
if n in minnodes:
nextmin=n
break
if nextmin>0:
npeaks=2
#if rtype==3:
# return nmax
if allup==True and npeaks==0:
return -1
if diss==True and npeaks==0:
return -2
return npeaks
geodesic_reparam(self)
inherited
🔗
Reparameterize using Geodesic interpolation
Source code in pygsm/growing_string_methods/de_gsm.py
def geodesic_reparam(self):
'''
Reparameterize using Geodesic interpolation
'''
printcool(' Reparameterizing using Geodesic Interpolation.')
TSnode = self.TSnode
if self.climb or self.find:
a = GSM.geodesic_reparam( self.nodes[0:self.TSnode] )
b = GSM.geodesic_reparam( self.nodes[self.TSnode:] )
new_xyzs = np.vstack((a,b))
else:
new_xyzs = GSM.geodesic_reparam(self.nodes)
for i,xyz in enumerate(new_xyzs):
self.nodes[i].xyz = xyz
self.refresh_coordinates()
get_tangent(node1, node2, print_level = 1, **kwargs)
inherited
🔗
Get internal coordinate tangent between two nodes, assumes they have unique IDs
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def get_tangent(node1,node2,print_level=1,**kwargs):
'''
Get internal coordinate tangent between two nodes, assumes they have unique IDs
'''
if node2 is not None and node1.node_id!=node2.node_id:
print(" getting tangent from between %i %i pointing towards %i"%(node2.node_id,node1.node_id,node2.node_id))
assert node2!=None,'node n2 is None'
PMDiff = np.zeros(node2.num_primitives)
for k,prim in enumerate(node2.primitive_internal_coordinates):
PMDiff[k] = prim.calcDiff(node2.xyz,node1.xyz)
return np.reshape(PMDiff,(-1,1)),None
else:
print(" getting tangent from node ",node1.node_id)
driving_coords = kwargs.get('driving_coords',None)
assert driving_coords is not None, " Driving coord is None!"
c = Counter(elem[0] for elem in driving_coords)
nadds = c['ADD']
nbreaks = c['BREAK']
nangles = c['nangles']
ntorsions = c['ntorsions']
ictan = np.zeros((node1.num_primitives,1),dtype=float)
breakdq = 0.3
bdist=0.0
atoms = node1.atoms
xyz = node1.xyz.copy()
for i in driving_coords:
if "ADD" in i:
#order indices to avoid duplicate bonds
if i[1]<i[2]:
index = [i[1]-1, i[2]-1]
else:
index = [i[2]-1, i[1]-1]
bond = Distance(index[0],index[1])
prim_idx = node1.coord_obj.Prims.dof_index(bond)
if len(i)==3:
#TODO why not just use the covalent radii?
d0 = (atoms[index[0]].vdw_radius + atoms[index[1]].vdw_radius)/2.8
elif len(i)==4:
d0=i[3]
current_d = bond.value(xyz)
#TODO don't set tangent if value is too small
ictan[prim_idx] = -1*(d0-current_d)
#if nbreaks>0:
# ictan[prim_idx] *= 2
# => calc bdist <=
if current_d>d0:
bdist += np.dot(ictan[prim_idx],ictan[prim_idx])
if print_level>0:
print(" bond %s target (less than): %4.3f current d: %4.3f diff: %4.3f " % ((i[1],i[2]),d0,current_d,ictan[prim_idx]))
elif "BREAK" in i:
#order indices to avoid duplicate bonds
if i[1]<i[2]:
index = [i[1]-1, i[2]-1]
else:
index = [i[2]-1, i[1]-1]
bond = Distance(index[0],index[1])
prim_idx = node1.coord_obj.Prims.dof_index(bond)
if len(i)==3:
d0 = (atoms[index[0]].vdw_radius + atoms[index[1]].vdw_radius)
elif len(i)==4:
d0=i[3]
current_d = bond.value(xyz)
ictan[prim_idx] = -1*(d0-current_d)
# => calc bdist <=
if current_d<d0:
bdist += np.dot(ictan[prim_idx],ictan[prim_idx])
if print_level>0:
print(" bond %s target (greater than): %4.3f, current d: %4.3f diff: %4.3f " % ((i[1],i[2]),d0,current_d,ictan[prim_idx]))
elif "ANGLE" in i:
if i[1]<i[3]:
index = [i[1]-1, i[2]-1,i[3]-1]
else:
index = [i[3]-1, i[2]-1,i[1]-1]
angle = Angle(index[0],index[1],index[2])
prim_idx = node1.coord_obj.Prims.dof_index(angle)
anglet = i[4]
ang_value = angle.value(xyz)
ang_diff = anglet*np.pi/180. - ang_value
#print(" angle: %s is index %i " %(angle,ang_idx))
if print_level>0:
print((" anglev: %4.3f align to %4.3f diff(rad): %4.3f" %(ang_value,anglet,ang_diff)))
ictan[prim_idx] = -ang_diff
#TODO need to come up with an adist
#if abs(ang_diff)>0.1:
# bdist+=ictan[ICoord1.BObj.nbonds+ang_idx]*ictan[ICoord1.BObj.nbonds+ang_idx]
elif "TORSION" in i:
if i[1]<i[4]:
index = [i[1]-1,i[2]-1,i[3]-1,i[4]-1]
else:
index = [i[4]-1,i[3]-1,i[2]-1,i[1]-1]
torsion = Dihedral(index[0],index[1],index[2],index[3])
prim_idx = node1.coord_obj.Prims.dof_index(torsion)
tort = i[5]
torv = torsion.value(xyz)
tor_diff = tort - torv*180./np.pi
if tor_diff>180.:
tor_diff-=360.
elif tor_diff<-180.:
tor_diff+=360.
ictan[prim_idx] = -tor_diff*np.pi/180.
if tor_diff*np.pi/180.>0.1 or tor_diff*np.pi/180.<0.1:
bdist += np.dot(ictan[prim_idx],ictan[prim_idx])
if print_level>0:
print((" current torv: %4.3f align to %4.3f diff(deg): %4.3f" %(torv*180./np.pi,tort,tor_diff)))
elif "OOP" in i:
index = [i[1]-1,i[2]-1,i[3]-1,i[4]-1]
oop = OutOfPlane(index[0],index[1],index[2],index[3])
prim_idx = node1.coord_obj.Prims.dof_index(oop)
oopt = i[5]
oopv = oop.value(xyz)
oop_diff = oopt - oopv*180./np.pi
if oop_diff>180.:
oop_diff-=360.
elif oop_diff<-180.:
oop_diff+=360.
ictan[prim_idx] = -oop_diff*np.pi/180.
if oop_diff*np.pi/180.>0.1 or oop_diff*np.pi/180.<0.1:
bdist += np.dot(ictan[prim_idx],ictan[prim_idx])
if print_level>0:
print((" current oopv: %4.3f align to %4.3f diff(deg): %4.3f" %(oopv*180./np.pi,oopt,oop_diff)))
bdist = np.sqrt(bdist)
if np.all(ictan==0.0):
raise RuntimeError(" All elements are zero")
return ictan,bdist
get_tangents(nodes, n0 = 0, print_level = 0)
inherited
🔗
Get the normalized internal coordinate tangents and magnitudes between all nodes
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def get_tangents(nodes,n0=0,print_level=0):
'''
Get the normalized internal coordinate tangents and magnitudes between all nodes
'''
nnodes = len(nodes)
dqmaga = [0.]*nnodes
dqa = np.zeros((nnodes+1,nnodes))
ictan = [[]]*nnodes
for n in range(n0+1,nnodes):
#print "getting tangent between %i %i" % (n,n-1)
assert nodes[n]!=None,"n is bad"
assert nodes[n-1]!=None,"n-1 is bad"
ictan[n],_ = GSM.get_tangent(nodes[n-1],nodes[n])
dqmaga[n] = 0.
#ictan0= np.copy(ictan[n])
dqmaga[n] = np.linalg.norm(ictan[n])
ictan[n] /= dqmaga[n]
# NOTE:
# vanilla GSM has a strange metric for distance
# no longer following 7/1/2020
#constraint = self.newic.constraints[:,0]
# just a fancy way to get the normalized tangent vector
#prim_constraint = block_matrix.dot(Vecs,constraint)
#for prim in self.newic.primitive_internal_coordinates:
# if type(prim) is Distance:
# index = self.newic.coord_obj.Prims.dof_index(prim)
# prim_constraint[index] *= 2.5
#dqmaga[n] = float(np.dot(prim_constraint.T,ictan0))
#dqmaga[n] = float(np.sqrt(dqmaga[n]))
if dqmaga[n]<0.:
raise RuntimeError
# TEMPORORARY parallel idea
#ictan = [0.]
#ictan += [ Process(target=get_tangent,args=(n,)) for n in range(n0+1,self.nnodes)]
#dqmaga = [ Process(target=get_dqmag,args=(n,ictan[n])) for n in range(n0+1,self.nnodes)]
if print_level>1:
print('------------printing ictan[:]-------------')
for n in range(n0+1,nnodes):
print("ictan[%i]" %n)
print(ictan[n].T)
if print_level>0:
print('------------printing dqmaga---------------')
for n in range(n0+1,nnodes):
print(" {:5.4}".format(dqmaga[n]),end='')
if (n)%5==0:
print()
print()
return ictan,dqmaga
get_tangents_growing(self, print_level = 1)
inherited
🔗
Finds the tangents during the growth phase. Tangents referenced to left or right during growing phase. Also updates coordinates Not a static method beause no one should ever call this outside of GSM
Source code in pygsm/growing_string_methods/de_gsm.py
def get_tangents_growing(self,print_level=1):
"""
Finds the tangents during the growth phase.
Tangents referenced to left or right during growing phase.
Also updates coordinates
Not a static method beause no one should ever call this outside of GSM
"""
ncurrent,nlist = self.make_difference_node_list()
dqmaga = [0.]*self.nnodes
ictan = [[]]*self.nnodes
if self.print_level>1:
print("ncurrent, nlist")
print(ncurrent)
print(nlist)
for n in range(ncurrent):
print(" ictan[{}]".format(nlist[2*n]))
ictan0,_ = self.get_tangent(
node1=self.nodes[nlist[2*n]],
node2=self.nodes[nlist[2*n+1]],
driving_coords=self.driving_coords,
)
if self.print_level>1:
print("forming space for", nlist[2*n+1])
if self.print_level>1:
print("forming tangent for ",nlist[2*n])
if (ictan0[:]==0.).all():
print(" ICTAN IS ZERO!")
print(nlist[2*n])
print(nlist[2*n+1])
raise RuntimeError
#normalize ictan
norm = np.linalg.norm(ictan0)
ictan[nlist[2*n]] = ictan0/norm
# NOTE regular GSM does something weird here
#Vecs = self.nodes[nlist[2*n]].update_coordinate_basis(constraints=self.ictan[nlist[2*n]])
#constraint = self.nodes[nlist[2*n]].constraints
#prim_constraint = block_matrix.dot(Vecs,constraint)
# but this is not followed here anymore 7/1/2020
#dqmaga[nlist[2*n]] = np.dot(prim_constraint.T,ictan0)
#dqmaga[nlist[2*n]] = float(np.sqrt(abs(dqmaga[nlist[2*n]])))
#tmp_dqmaga = np.dot(prim_constraint.T,ictan0)
#tmp_dqmaga = np.sqrt(tmp_dqmaga)
dqmaga[nlist[2*n]] = norm
if print_level>0:
print('------------printing dqmaga---------------')
for n in range(self.nnodes):
print(" {:5.3}".format(dqmaga[n]), end=' ')
if (n+1)%5==0:
print()
print()
if print_level>1:
for n in range(ncurrent):
print("dqmag[%i] =%1.2f" %(nlist[2*n],self.dqmaga[nlist[2*n]]))
print("printing ictan[%i]" %nlist[2*n])
print(self.ictan[nlist[2*n]].T)
for i,tan in enumerate(ictan):
if np.all(tan==0.0):
print("tan %i of the tangents is 0" %i)
raise RuntimeError
return ictan,dqmaga
get_three_way_tangents(nodes, energies, find = True, n0 = 0, print_level = 0)
inherited
🔗
Calculates internal coordinate tangent with a three-way tangent at TS node
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def get_three_way_tangents(nodes,energies,find=True,n0=0,print_level=0):
'''
Calculates internal coordinate tangent with a three-way tangent at TS node
'''
nnodes = len(nodes)
ictan = [[]]*nnodes
dqmaga = [0.]*nnodes
#TSnode = np.argmax(energies[1:nnodes-1])+1
TSnode = np.argmax(energies) # allow for the possibility of TS node to be endpoints?
last_node_max = (TSnode==nnodes-1)
first_node_max = (TSnode==0)
if first_node_max or last_node_max:
print("*********** This will cause a range error in the following for loop *********")
print("** Setting the middle of the string to be TS node to get proper directions **")
TSnode = nnodes//2
for n in range(n0+1,nnodes):
do3 = False
#if not find:
# if energies[n+1] > energies[n] and energies[n] > energies[n-1]:
# intic_n = n
# newic_n = n+1
# elif energies[n-1] > energies[n] and energies[n] > energies[n+1]:
# #intic_n = n-1
# #newic_n = n
# intic_n = n
# newic_n = n-1
# else:
# do3 = True
# newic_n = n
# intic_n = n+1
# int2ic_n = n-1
#else:
if n < TSnode:
# The order is very important here
intic_n = n
newic_n = n+1
elif n> TSnode:
# The order is very important here
intic_n = n
newic_n = n-1
else:
do3 = True
newic_n = n
intic_n = n+1
int2ic_n = n-1
if do3:
if first_node_max or last_node_max:
t1,_ = GSM.get_tangent(nodes[intic_n],nodes[newic_n])
t2,_ = GSM.get_tangent(nodes[newic_n],nodes[int2ic_n])
print(" done 3 way tangent")
ictan0 = t1 + t2
else:
f1 = 0.
dE1 = abs(energies[n+1]-energies[n])
dE2 = abs(energies[n] - energies[n-1])
dEmax = max(dE1,dE2)
dEmin = min(dE1,dE2)
if energies[n+1]>energies[n-1]:
f1 = dEmax/(dEmax+dEmin+0.00000001)
else:
f1 = 1 - dEmax/(dEmax+dEmin+0.00000001)
print(' 3 way tangent ({}): f1:{:3.2}'.format(n,f1))
t1,_ = GSM.get_tangent(nodes[intic_n],nodes[newic_n])
t2,_ = GSM.get_tangent(nodes[newic_n],nodes[int2ic_n])
print(" done 3 way tangent")
ictan0 = f1*t1 +(1.-f1)*t2
else:
ictan0,_ = GSM.get_tangent(nodes[newic_n],nodes[intic_n])
ictan[n] = ictan0/np.linalg.norm(ictan0)
dqmaga[n]=np.linalg.norm(ictan0)
return ictan,dqmaga
go_gsm(self, max_iters = 50, opt_steps = 3, rtype = 2)
🔗
rtype=2 Find and Climb TS, 1 Climb with no exact find, 0 turning of climbing image and TS search
Source code in pygsm/growing_string_methods/de_gsm.py
def go_gsm(self,max_iters=50,opt_steps=3,rtype=2):
"""
rtype=2 Find and Climb TS,
1 Climb with no exact find,
0 turning of climbing image and TS search
"""
self.set_V0()
if not self.isRestarted:
if self.growth_direction==0:
self.add_GSM_nodes(2)
elif self.growth_direction==1:
self.add_GSM_nodeR(1)
elif self.growth_direction==2:
self.add_GSM_nodeP(1)
# Grow String
self.grow_string(max_iters=max_iters,max_opt_steps=opt_steps)
nifty.printcool("Done Growing the String!!!")
self.done_growing = True
#nifty.printcool("initial ic_reparam")
self.reparameterize()
write_molden_geoms('grown_string_{:03}.xyz'.format(self.ID),self.geometries,self.energies,self.gradrmss,self.dEs)
else:
self.ictan,self.dqmaga = self.get_tangents(self.nodes)
self.refresh_coordinates()
if self.tscontinue:
self.optimize_string(max_iter=max_iters,opt_steps=opt_steps,rtype=rtype)
else:
print("Exiting early")
self.end_early=True
print("Finished GSM!")
return self.nnodes,self.energies
grow_nodes(self)
🔗
Grow nodes
Source code in pygsm/growing_string_methods/de_gsm.py
def grow_nodes(self):
'''
Grow nodes
'''
if self.nodes[self.nR-1].gradrms < self.gaddmax and self.growth_direction!=2:
if self.nodes[self.nR] == None:
self.add_GSM_nodeR()
print(" getting energy for node %d: %5.4f" %(self.nR-1,self.nodes[self.nR-1].energy - self.nodes[0].V0))
if self.nodes[self.nnodes-self.nP].gradrms < self.gaddmax and self.growth_direction!=1:
if self.nodes[-self.nP-1] == None:
self.add_GSM_nodeP()
print(" getting energy for node %d: %5.4f" %(self.nnodes-self.nP,self.nodes[-self.nP].energy - self.nodes[0].V0))
return
grow_string(self, max_iters = 30, max_opt_steps = 3, nconstraints = 1)
inherited
🔗
Grow the string
Parameters🔗
max_iter : int Maximum number of GSM iterations nconstraints : int optsteps : int Maximum number of optimization steps per node of string
Source code in pygsm/growing_string_methods/de_gsm.py
def grow_string(self,max_iters=30,max_opt_steps=3,nconstraints=1):
'''
Grow the string
Parameters
----------
max_iter : int
Maximum number of GSM iterations
nconstraints : int
optsteps : int
Maximum number of optimization steps per node of string
'''
printcool("In growth_iters")
ncurrent,nlist = self.make_difference_node_list()
self.ictan,self.dqmaga = self.get_tangents_growing()
self.refresh_coordinates()
self.set_active(self.nR-1, self.nnodes-self.nP)
isGrown=False
iteration=0
while not isGrown:
if iteration>max_iters:
raise RuntimeError
printcool("Starting growth iteration %i" % iteration)
self.optimize_iteration(max_opt_steps)
totalgrad,gradrms,sum_gradrms = self.calc_optimization_metrics(self.nodes)
write_molden_geoms('scratch/growth_iters_{:03}_{:03}.xyz'.format(self.ID,iteration),self.geometries,self.energies,self.gradrmss,self.dEs)
print(" gopt_iter: {:2} totalgrad: {:4.3} gradrms: {:5.4} max E: {:5.4}\n".format(iteration,float(totalgrad),float(gradrms),float(self.emax)))
try:
self.grow_nodes()
except Exception as error:
print("can't add anymore nodes, bdist too small")
if self.__class__.__name__=="SE_GSM": # or self.__class__.__name__=="SE_Cross":
# Don't do SE_cross because that already does optimization later
if self.nodes[self.nR-1].PES.lot.do_coupling:
opt_type='MECI'
else:
opt_type='UNCONSTRAINED'
print(" optimizing last node")
self.optimizer[self.nR-1].conv_grms = self.CONV_TOL
print(self.optimizer[self.nR-1].conv_grms)
path=os.path.join(os.getcwd(),'scratch/{:03d}/{}'.format(self.ID,self.nR-1))
self.optimizer[self.nR-1].optimize(
molecule=self.nodes[self.nR-1],
refE=self.nodes[0].V0,
opt_steps=50,
opt_type=opt_type,
path=path,
)
elif self.__class__.__name__=="SE_Cross":
print(" Will do extra optimization of this node in SE-Cross")
else:
raise RuntimeError
self.set_active(self.nR-1, self.nnodes-self.nP)
self.ic_reparam_g()
self.ictan,self.dqmaga = self.get_tangents_growing()
self.refresh_coordinates()
iteration+=1
isGrown = self.check_if_grown()
# create newic object
print(" creating newic molecule--used for ic_reparam")
self.newic = Molecule.copy_from_options(self.nodes[0])
# TODO should something be done for growthdirection 2?
if self.growth_direction==1:
print("Setting LOT of last node")
self.nodes[-1] = Molecule.copy_from_options(
MoleculeA = self.nodes[-2],
xyz = self.nodes[-1].xyz,
new_node_id = self.nnodes-1
)
return
has_intermediate(self)
inherited
🔗
Check string for intermediates
Source code in pygsm/growing_string_methods/de_gsm.py
def has_intermediate(self):
'''
Check string for intermediates
'''
energies = self.energies
potential_min = []
for i in range(1, (len(energies) - 1)):
rnoise = 0
pnoise = 0
a = 1
b = 1
while (energies[i-a] >= energies[i]):
if (energies[i-a] - energies[i]) > rnoise:
rnoise = energies[i-a] - energies[i]
if rnoise > self.noise:
break
if (i-a) == 0:
break
a += 1
while (energies[i+b] >= energies[i]):
if (energies[i+b] - energies[i]) > pnoise:
pnoise = energies[i+b] - energies[i]
if pnoise > self.noise:
break
if (i+b) == len(energies) - 1:
break
b += 1
if ((rnoise > self.noise) and (pnoise > self.noise)):
print('Potential minimum at image %s', str(i))
potential_min.append(i)
return len(potential_min)>0
ic_reparam(nodes, energies, climbing = False, ic_reparam_steps = 8, print_level = 1)
inherited
🔗
Reparameterizes the string using Delocalizedin internal coordinatesusing three-way tangents at the TS node Only pushes nodes outwards during reparameterization because otherwise too many things change.
Be careful, however, if the path is allup or alldown then this can cause
Parameters🔗
nodes : list of molecule objects energies : list of energies in kcal/mol ic_reparam_steps : int max number of reparameterization steps print_level : int verbosity
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def ic_reparam(nodes,energies,climbing=False,ic_reparam_steps=8,print_level=1):
'''
Reparameterizes the string using Delocalizedin internal coordinatesusing three-way tangents at the TS node
Only pushes nodes outwards during reparameterization because otherwise too many things change.
Be careful, however, if the path is allup or alldown then this can cause
Parameters
----------
nodes : list of molecule objects
energies : list of energies in kcal/mol
ic_reparam_steps : int max number of reparameterization steps
print_level : int verbosity
'''
nifty.printcool("reparametrizing string nodes")
nnodes=len(nodes)
rpart = np.zeros(nnodes)
for n in range(1,nnodes):
rpart[n] = 1./(nnodes-1)
deltadqs = np.zeros(nnodes)
TSnode = np.argmax(energies)
MAXRE = 0.5
ictan,orig_dqmaga = GSM.get_three_way_tangents(nodes,energies)
if (TSnode==nnodes-1) or (TSnode==0):
print("********** This will cause a range error in the following for loop *********")
print("** Setting the middle of the string to be TS node to get proper directions **")
TSnode = nnodes//2
if climbing:
raise RuntimeError(" This shouldn't happen")
ideal_progress_gained = np.zeros(nnodes)
if climbing:
for n in range(1,TSnode):
ideal_progress_gained[n] = 1./(TSnode-1)
for n in range(TSnode+1,nnodes):
ideal_progress_gained[n] = 1./(nnodes-TSnode-1)
ideal_progress_gained[TSnode]=0.
else:
for n in range(1,nnodes):
ideal_progress_gained[n] = 1./(nnodes-1)
for i in range(ic_reparam_steps):
ictan,dqmaga = GSM.get_three_way_tangents(nodes,energies)
totaldqmag = np.sum(dqmaga[1:nnodes])
if climbing:
progress=np.zeros_like(dqmaga)
progress_gained=np.zeros_like(dqmaga)
h1dqmag = np.sum(dqmaga[1:TSnode])
h2dqmag = np.sum(dqmaga[TSnode+1:nnodes])
if print_level>0:
print(" h1dqmag, h2dqmag: %3.2f %3.2f" % (h1dqmag,h2dqmag))
progress_gained[1:TSnode] = dqmaga[1:TSnode]/h1dqmag
progress_gained[TSnode+1:] = dqmaga[TSnode+1:]/h2dqmag
progress[1:TSnode] = np.cumsum(progress_gained[1:TSnode])
progress[TSnode+1:] = np.cumsum(progress_gained[TSnode+1:])
else:
progress = np.cumsum(dqmaga)/totaldqmag
progress_gained = dqmaga/totaldqmag
if i==0:
orig_progress_gained = copy(progress_gained)
if climbing:
difference=np.zeros_like(dqmaga)
for n in range(1,TSnode):
difference[n] = ideal_progress_gained[n] - progress_gained[n]
deltadqs[n] = difference[n]*h1dqmag
for n in range(TSnode+1,nnodes):
difference[n] = ideal_progress_gained[n] - progress_gained[n]
deltadqs[n] = difference[n]*h2dqmag
else:
difference = ideal_progress_gained - progress_gained
deltadqs = difference*totaldqmag
if print_level>1:
print('dqmaga')
print(dqmaga[1:],len(dqmaga))
print(" progress")
print(progress)
print(" progress gained per step")
print(progress_gained)
print('difference')
print(difference)
print("move each node this much")
print(deltadqs)
# Move left wing outwards only
for n in range(1,TSnode+1):
# Move along to previous tan
if deltadqs[n]<0:
print(f" Moving node {n} along tan[{n}] this much {deltadqs[n]}")
if climbing and n==TSnode:
assert deltadqs[n]==0, "This is wrong."
nodes[n].update_coordinate_basis(ictan[n])
constraint = nodes[n].constraints[:,0]
dq = deltadqs[n]*constraint
nodes[n].update_xyz(dq,verbose=(print_level>1))
# Move right wing outwards only
for n in range(TSnode,nnodes-1):
# Move along to next tan
if deltadqs[n+1]<0:
if climbing and n==TSnode:
print('skipping TS')
continue
print(f" Moving node {n} along tan[{n+1}] this much {deltadqs[n+1]}")
nodes[n].update_coordinate_basis(ictan[n+1])
constraint = nodes[n].constraints[:,0]
dq = -deltadqs[n+1]*constraint
nodes[n].update_xyz(dq,verbose=(print_level>1))
disprms = np.linalg.norm(deltadqs)/np.sqrt(nnodes-1)
if disprms < 0.02:
break
if print_level>0:
for n in range(1,nnodes-1):
print(" disp[{}]: {:1.3}".format(n,deltadqs[n]), end=' ')
print()
print(" disprms: {:1.3}\n".format(disprms))
ictan,dqmaga = GSM.get_three_way_tangents(nodes,energies)
print('\n')
if print_level>0:
print(" ideal progress gained per step",end=' ')
for n in range(1,nnodes):
print(" step [{}]: {:1.3}".format(n,ideal_progress_gained[n]), end=' ')
print()
print(" original path progress",end=' ' )
for n in range(1,nnodes):
print(" step [{}]: {:1.3}".format(n,orig_progress_gained[n]), end=' ')
print()
print(" reparameterized path progress",end=' ')
for n in range(1,nnodes):
print(" step [{}]: {:1.3}".format(n,progress_gained[n]), end=' ')
print()
print(" spacings (end ic_reparam, steps: {}/{}):".format(i+1,ic_reparam_steps), end=' ')
for n in range(nnodes):
print(" {:1.2}".format(dqmaga[n]), end=' ')
print("\n disprms: {:1.3}".format(disprms))
return
ic_reparam_g(self, ic_reparam_steps = 4, n0 = 0, reparam_interior = True)
inherited
🔗
Reparameterize during growth phase
Source code in pygsm/growing_string_methods/de_gsm.py
def ic_reparam_g(self,ic_reparam_steps=4,n0=0,reparam_interior=True): #see line 3863 of gstring.cpp
"""
Reparameterize during growth phase
"""
printcool("Reparamerizing string nodes")
#close_dist_fix(0) #done here in GString line 3427.
rpmove = np.zeros(self.nnodes)
rpart = np.zeros(self.nnodes)
dqavg = 0.0
disprms = 0.0
h1dqmag = 0.0
h2dqmag = 0.0
dE = np.zeros(self.nnodes)
edist = np.zeros(self.nnodes)
emax = -1000 # And this?
if self.current_nnodes==self.nnodes:
self.ic_reparam(4)
return
for i in range(ic_reparam_steps):
self.ictan,self.dqmaga = self.get_tangents_growing()
self.refresh_coordinates()
totaldqmag = np.sum(self.dqmaga[n0:self.nR-1])+np.sum(self.dqmaga[self.nnodes-self.nP+1:self.nnodes])
if self.print_level>0:
if i==0:
print(" totaldqmag (without inner): {:1.2}\n".format(totaldqmag))
print(" printing spacings dqmaga: ")
for n in range(self.nnodes):
print(" {:2.3}".format(self.dqmaga[n]), end=' ')
if (n+1)%5==0:
print()
print()
if i == 0:
if self.current_nnodes!=self.nnodes:
rpart = np.zeros(self.nnodes)
for n in range(n0+1,self.nR):
rpart[n] = 1.0/(self.current_nnodes-2)
for n in range(self.nnodes-self.nP,self.nnodes-1):
rpart[n] = 1.0/(self.current_nnodes-2)
else:
for n in range(n0+1,self.nnodes):
rpart[n] = 1./(self.nnodes-1)
if self.print_level>0:
if i==0:
print(" rpart: ")
for n in range(1,self.nnodes-1):
print(" {:1.2}".format(rpart[n]), end=' ')
if (n)%5==0:
print()
print()
nR0 = self.nR
nP0 = self.nP
# TODO CRA 3/2019 why is this here?
if not reparam_interior:
if self.nnodes-self.current_nnodes > 2:
nR0 -= 1
nP0 -= 1
deltadq = 0.0
for n in range(n0+1,nR0):
deltadq = self.dqmaga[n-1] - totaldqmag*rpart[n]
rpmove[n] = -deltadq
for n in range(self.nnodes-nP0,self.nnodes-1):
deltadq = self.dqmaga[n+1] - totaldqmag*rpart[n]
rpmove[n] = -deltadq
MAXRE = 1.1
for n in range(n0+1,self.nnodes-1):
if abs(rpmove[n]) > MAXRE:
rpmove[n] = float(np.sign(rpmove[n])*MAXRE)
disprms = float(np.linalg.norm(rpmove[n0+1:self.nnodes-1]))
lastdispr = disprms
if self.print_level>0:
for n in range(n0+1,self.nnodes-1):
print(" disp[{}]: {:1.2f}".format(n,rpmove[n]), end=' ')
if (n)%5==0:
print()
print()
print(" disprms: {:1.3}\n".format(disprms))
if disprms < 1e-2:
break
ncurrent,nlist = self.make_difference_node_list()
param_list=[]
for n in range(ncurrent):
if nlist[2*n+1] not in param_list:
if rpmove[nlist[2*n+1]]>0:
# Using tangent pointing inner?
print('Moving {} along ictan[{}]'.format(nlist[2*n+1],nlist[2*n+1]))
self.nodes[nlist[2*n+1]].update_coordinate_basis(constraints=self.ictan[nlist[2*n+1]])
constraint = self.nodes[nlist[2*n+1]].constraints[:,0]
dq0 = rpmove[nlist[2*n+1]]*constraint
self.nodes[nlist[2*n+1]].update_xyz(dq0,verbose=True)
param_list.append(nlist[2*n+1])
else:
# Using tangent point outer
print('Moving {} along ictan[{}]'.format(nlist[2*n+1],nlist[2*n]))
self.nodes[nlist[2*n+1]].update_coordinate_basis(constraints=self.ictan[nlist[2*n]])
constraint = self.nodes[nlist[2*n+1]].constraints[:,0]
dq0 = rpmove[nlist[2*n+1]]*constraint
self.nodes[nlist[2*n+1]].update_xyz(dq0,verbose=True)
param_list.append(nlist[2*n+1])
print(" spacings (end ic_reparam, steps: {}/{}):".format(i+1,ic_reparam_steps), end=' ')
for n in range(self.nnodes):
print(" {:1.2}".format(self.dqmaga[n]), end=' ')
print(" disprms: {:1.3}".format(disprms))
#TODO old GSM does this here
#Failed = check_array(self.nnodes,self.dqmaga)
#If failed, do exit 1
interpolate(start_node, end_node, num_interp)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def interpolate(start_node,end_node,num_interp):
'''
'''
nifty.printcool(" interpolate")
num_nodes = num_interp + 2
nodes = [None]*(num_nodes)
nodes[0] = start_node
nodes[-1] = end_node
sign=1
nR = 1
nP = 1
nn = nR + nP
for n in range(num_interp):
if num_nodes - nn > 1:
stepsize = 1./float(num_nodes - nn)
else:
stepsize = 0.5
if sign == 1:
iR = nR-1
iP = num_nodes - nP
iN = nR
nodes[nR] = GSM.add_node(nodes[iR],nodes[iP],stepsize,iN)
if nodes[nR] == None:
raise RuntimeError
#print(" Energy of node {} is {:5.4}".format(nR,nodes[nR].energy-E0))
nR +=1
nn += 1
else:
n1 = num_nodes - nP
n2 = n1 -1
n3 = nR -1
nodes[n2] = GSM.add_node(nodes[n1],nodes[n3],stepsize,n2)
if nodes[n2] == None:
raise RuntimeError
#print(" Energy of node {} is {:5.4}".format(nR,nodes[nR].energy-E0))
nP +=1
nn += 1
sign *= -1
return nodes
interpolate_xyz(nodeR, nodeP, stepsize)
inherited
🔗
Interpolate between two nodes
Source code in pygsm/growing_string_methods/de_gsm.py
@staticmethod
def interpolate_xyz(nodeR,nodeP,stepsize):
'''
Interpolate between two nodes
'''
ictan,_ = get_tangent(nodeR,nodeP)
Vecs = nodeR.update_coordinate_basis(constraints=ictan)
constraint = nodeR.constraints[:,0]
prim_constraint = block_matrix.dot(Vecs,constraint)
dqmag = np.dot(prim_constraint.T,ictan)
print(" dqmag: %1.3f"%dqmag)
#sign=-1
sign=1.
dqmag *= (sign*stepsize)
print(" scaled dqmag: %1.3f"%dqmag)
dq0 = dqmag*constraint
old_xyz = nodeR.xyz.copy()
new_xyz = nodeR.coord_obj.newCartesian(old_xyz,dq0)
return new_xyz
is_converged(self, totalgrad, fp, rtype, ts_cgradq)
inherited
🔗
Check if optimization is converged
Source code in pygsm/growing_string_methods/de_gsm.py
def is_converged(self,totalgrad,fp,rtype,ts_cgradq):
'''
Check if optimization is converged
'''
TS_conv = self.options['CONV_TOL']
#if self.find and self.optimizer[self.TSnode].nneg>1:
# print(" reducing TS convergence because nneg>1")
# TS_conv = self.options['CONV_TOL']/2.
self.optimizer[self.TSnode].conv_grms = TS_conv
#print(" Number of imaginary frequencies %i" % self.optimizer[self.TSnode].nneg)
if (rtype == 2 and self.find):
return (self.nodes[self.TSnode].gradrms<TS_conv and self.dE_iter < self.optimizer[self.TSnode].conv_Ediff) or (totalgrad<0.1 and self.nodes[self.TSnode].gradrms<2.5*TS_conv and self.dE_iter<0.02 and self.optimizer[self.TSnode].nneg <2) #TODO extra crit here
elif rtype==1 and self.climb:
return self.nodes[self.TSnode].gradrms<TS_conv and abs(ts_cgradq) < self.options['CONV_TOL'] and self.dE_iter < 0.2
elif not self.climber and not self.finder:
print(" CONV_TOL=%.4f" %self.CONV_TOL)
return all([self.optimizer[n].converged for n in range(self.nnodes)])
# => Check if intermediate exists
if self.has_intermediate():
self.endearly=True #bools
self.tscontinue=False
return True
return False
make_difference_node_list(self)
🔗
Returns ncurrent and a list of indices that can be iterated over to produce tangents for the string pathway.
Source code in pygsm/growing_string_methods/de_gsm.py
def make_difference_node_list(self):
'''
Returns ncurrent and a list of indices that can be iterated over to produce
tangents for the string pathway.
'''
#TODO: THis can probably be done more succinctly using a list of tuples
ncurrent = 0
nlist = [0]*(2*self.nnodes)
for n in range(self.nR-1):
nlist[2*ncurrent] = n
nlist[2*ncurrent+1] = n+1
ncurrent += 1
for n in range(self.nnodes-self.nP+1,self.nnodes):
nlist[2*ncurrent] = n
nlist[2*ncurrent+1] = n-1
ncurrent += 1
nlist[2*ncurrent] = self.nR -1
nlist[2*ncurrent+1] = self.nnodes - self.nP
if False:
nlist[2*ncurrent+1] = self.nR - 2 #for isMAP_SE
#TODO is this actually used?
#if self.nR == 0: nlist[2*ncurrent] += 1
#if self.nP == 0: nlist[2*ncurrent+1] -= 1
ncurrent += 1
nlist[2*ncurrent] = self.nnodes -self.nP
nlist[2*ncurrent+1] = self.nR-1
##TODO is this actually used?
#if self.nR == 0: nlist[2*ncurrent+1] += 1
#if self.nP == 0: nlist[2*ncurrent] -= 1
ncurrent += 1
return ncurrent,nlist
modify_TS_Hess(self)
inherited
🔗
Modifies Hessian using RP direction
Source code in pygsm/growing_string_methods/de_gsm.py
def modify_TS_Hess(self):
''' Modifies Hessian using RP direction'''
print("modifying %i Hessian with RP" % self.TSnode)
TSnode = self.TSnode
# a variable to determine how many time since last modify
self.hess_counter = 0
self.TS_E_0 = self.energies[self.TSnode]
E0 = self.energies[TSnode]/GSM.units.KCAL_MOL_PER_AU
Em1 = self.energies[TSnode-1]/GSM.units.KCAL_MOL_PER_AU
if self.TSnode+1<self.nnodes:
Ep1 = self.energies[TSnode+1]/GSM.units.KCAL_MOL_PER_AU
else:
Ep1 = Em1
# Update TS node coord basis
Vecs = self.nodes[TSnode].update_coordinate_basis(constraints=None)
# get constrained coord basis
self.newic.xyz = self.nodes[TSnode].xyz.copy()
const_vec = self.newic.update_coordinate_basis(constraints=self.ictan[TSnode])
q0 = self.newic.coordinates[0]
constraint = self.newic.constraints[:,0]
# this should just give back ictan[TSnode]?
tan0 = block_matrix.dot(const_vec,constraint)
# get qm1 (don't update basis)
self.newic.xyz = self.nodes[TSnode-1].xyz.copy()
qm1 = self.newic.coordinates[0]
if TSnode+1<self.nnodes:
# get qp1 (don't update basis)
self.newic.xyz = self.nodes[TSnode+1].xyz.copy()
qp1 = self.newic.coordinates[0]
else:
qp1 = qm1
print(" TS Hess init'd w/ existing Hintp")
# Go to non-constrained basis
self.newic.xyz = self.nodes[TSnode].xyz.copy()
self.newic.coord_basis = Vecs
self.newic.Primitive_Hessian = self.nodes[TSnode].Primitive_Hessian.copy()
self.newic.form_Hessian_in_basis()
tan = block_matrix.dot(block_matrix.transpose(Vecs),tan0) # (nicd,1
Ht = np.dot(self.newic.Hessian,tan) # (nicd,nicd)(nicd,1) = nicd,1
tHt = np.dot(tan.T,Ht)
a = abs(q0-qm1)
b = abs(qp1-q0)
c = 2*(Em1/a/(a+b) - E0/a/b + Ep1/b/(a+b))
print(" tHt %1.3f a: %1.1f b: %1.1f c: %1.3f" % (tHt,a[0],b[0],c[0]))
ttt = np.outer(tan,tan)
# Hint before
#with np.printoptions(threshold=np.inf):
# print self.newic.Hessian
#eig,tmph = np.linalg.eigh(self.newic.Hessian)
#print "initial eigenvalues"
#print eig
# Finalize Hessian
self.newic.Hessian += (c-tHt)*ttt
self.nodes[TSnode].Hessian = self.newic.Hessian.copy()
# Hint after
#with np.printoptions(threshold=np.inf):
# print self.nodes[TSnode].Hessian
#print "shape of Hessian is %s" % (np.shape(self.nodes[TSnode].Hessian),)
self.nodes[TSnode].newHess = 5
if False:
print("newHess of node %i %i" % (TSnode,self.nodes[TSnode].newHess))
eigen,tmph = np.linalg.eigh(self.nodes[TSnode].Hessian) #nicd,nicd
print("eigenvalues of new Hess")
print(eigen)
# reset pgradrms ?
mult_steps(self, n, opt_steps)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def mult_steps(self,n,opt_steps):
exsteps=1
tsnode = int(self.TSnode)
if (self.find) and self.energies[n] > self.energies[self.TSnode]*0.9 and n!=tsnode: #
exsteps=2
print(" multiplying steps for node %i by %i" % (n,exsteps))
self.optimizer[n].conv_grms = self.options['CONV_TOL'] # TODO this is not perfect here
self.optimizer[n].conv_gmax = self.options['CONV_gmax']
self.optimizer[n].conv_Ediff = self.options['CONV_Ediff']
#if (self.find or (self.climb and self.energies[tsnode]>self.energies[tsnode-1]+5 and self.energies[tsnode]>self.energies[tsnode+1]+5.)) and n==tsnode: #or self.climb
# exsteps=2
# print(" multiplying steps for node %i by %i" % (n,exsteps))
#elif not (self.find and self.climb) and self.energies[tsnode] > 1.75*self.energies[tsnode-1] and self.energies[tsnode] > 1.75*self.energies[tsnode+1] and self.done_growing and n==tsnode: #or self.climb
# exsteps=2
# print(" multiplying steps for node %i by %i" % (n,exsteps))
return exsteps*opt_steps
optimize_iteration(self, opt_steps)
inherited
🔗
Optimize string iteration
Source code in pygsm/growing_string_methods/de_gsm.py
def optimize_iteration(self,opt_steps):
'''
Optimize string iteration
'''
refE=self.nodes[0].energy
if self.use_multiprocessing:
cpus = mp.cpu_count()/self.nodes[0].PES.lot.nproc
print(" Parallelizing over {} processes".format(cpus))
out_queue = mp.Queue()
workers = [ mp.Process(target=mod, args=(self.nodes[n],self.optimizer[n],self.ictan[n],self.mult_steps(n,opt_steps),self.set_opt_type(n),refE,n,s,gp_prim, out_queue) ) for n in range(self.nnodes) if self.nodes[n] and self.active[n] ]
for work in workers: work.start()
for work in workers: work.join()
res_lst = []
for j in range(len(workers)):
res_lst.append(out_queue.get())
else:
for n in range(self.nnodes):
if self.nodes[n] and self.active[n]:
print()
path=os.path.join(os.getcwd(),'scratch/{:03d}/{}'.format(self.ID,n))
printcool("Optimizing node {}".format(n))
opt_type = self.set_opt_type(n)
osteps = self.mult_steps(n,opt_steps)
self.optimizer[n].optimize(
molecule=self.nodes[n],
refE=refE,
opt_type=opt_type,
opt_steps=osteps,
ictan=self.ictan[n],
xyzframerate=1,
path=path,
)
if self.__class__.__name__=="SE-GSM" and self.done_growing:
fp = self.find_peaks(2)
if self.energies[self.nnodes-1]>self.energies[self.nnodes-2] and fp>0 and self.nodes[self.nnodes-1].gradrms>self.CONV_TOL:
path=os.path.join(os.getcwd(),'scratch/{:03d}/{}'.format(self.ID,self.nnodes-1))
self.optimizer[self.nnodes-1].optimize(
molecule=self.nodes[self.nnodes-1],
refE=refE,
opt_type='UNCONSTRAINED',
opt_steps=osteps,
ictan=None,
path=path
)
optimize_string(self, max_iter = 30, nconstraints = 1, opt_steps = 1, rtype = 2)
inherited
🔗
Optimize the grown string until convergence
Parameters🔗
max_iter : int Maximum number of GSM iterations nconstraints : int optsteps : int Maximum number of optimization steps per node of string rtype : int An option to change how GSM optimizes
Source code in pygsm/growing_string_methods/de_gsm.py
def optimize_string(self,max_iter=30,nconstraints=1,opt_steps=1,rtype=2):
'''
Optimize the grown string until convergence
Parameters
----------
max_iter : int
Maximum number of GSM iterations
nconstraints : int
optsteps : int
Maximum number of optimization steps per node of string
rtype : int
An option to change how GSM optimizes
'''
printcool("In opt_iters")
self.nclimb=0
self.nhessreset=10 # are these used??? TODO
self.hessrcount=0 # are these used?! TODO
self.newclimbscale=2.
self.set_finder(rtype)
# set convergence for nodes
if (self.climber or self.finder):
factor = 2.5
else:
factor = 1.
for i in range(self.nnodes):
if self.nodes[i] !=None:
self.optimizer[i].conv_grms = self.CONV_TOL*factor
self.optimizer[i].conv_gmax = self.options['CONV_gmax']*factor
self.optimizer[i].conv_Ediff = self.options['CONV_Ediff']*factor
# enter loop
for oi in range(max_iter):
printcool("Starting opt iter %i" % oi)
if self.climb and not self.find: print(" CLIMBING")
elif self.find: print(" TS SEARCHING")
# stash previous TSnode
self.pTSnode = self.TSnode
self.emaxp = self.emax
# => Get all tangents 3-way <= #
self.ictan,self.dqmaga = self.get_three_way_tangents(self.nodes,self.energies)
# => do opt steps <= #
self.optimize_iteration(opt_steps)
print(" V_profile: ", end=' ')
self.print_energies()
#TODO resetting
#TODO special SSM criteria if first opt'd node is too high?
if (self.TSnode == self.nnodes-2 or self.TSnode == self.nnodes-3) and (self.climb or self.find):
printcool("WARNING\n: TS node shouldn't be second to last node for tangent reasons")
new_node = GSM.add_node(
self.nodes[self.TSnode],
self.nodes[self.TSnode+1],
stepsize=0.5,
node_id = self.TSnode+1,
)
new_node_list = [None]*(self.nnodes+1)
new_optimizers = [None]*(self.nnodes+1)
for n in range(0,self.TSnode+1):
new_node_list[n] = self.nodes[n]
new_optimizers[n] = self.optimizer[n]
new_node_list[self.TSnode+1] = new_node
new_optimizers[self.TSnode+1] = self.optimizer[0].__class__(self.optimizer[0].options.copy())
for n in range(self.TSnode+2,self.nnodes+1):
new_node_list[n] = Molecule.copy_from_options(MoleculeA = self.nodes[n-1], new_node_id = n)
new_optimizers[n] = self.optimizer[n-1]
self.nodes = new_node_list
self.optimizer = new_optimizers
self.nnodes = len(self.nodes)
print('New number of nodes %d' % self.nnodes)
self.active = [True] * self.nnodes
self.active[0] = False
self.active[self.nnodes-1] = False
if (self.TSnode == 1 or self.TSnode==2) and (self.climb or self.find):
printcool("WARNING\n: TS node shouldn't be first or second node for tangent reasons")
new_node = GSM.add_node(
self.nodes[self.TSnode-1],
self.nodes[self.TSnode],
stepsize=0.5,
node_id = self.TSnode-1,
)
new_node_list = [None]*(self.nnodes+1)
new_optimizers = [None]*(self.nnodes+1)
for n in range(0,self.TSnode-1):
new_node_list[n] = self.nodes[n]
new_optimizers[n] = self.optimizer[n]
new_node_list[self.TSnode-1] = new_node
new_optimizers[self.TSnode-1] = self.optimizer[0].__class__(self.optimizer[0].options.copy())
for n in range(self.TSnode,self.nnodes+1):
new_node_list[n] = Molecule.copy_from_options(MoleculeA = self.nodes[n-1], new_node_id = n)
new_optimizers[n] = self.optimizer[n-1]
self.nodes = new_node_list
self.optimizer = new_optimizers
self.nnodes = len(self.nodes)
print('New number of nodes %d' % self.nnodes)
self.active = [True] * self.nnodes
self.active[0] = False
self.active[self.nnodes-1] = False
# => find peaks <= #
fp = self.find_peaks(2)
ts_cgradq = 0.
if not self.find:
ts_cgradq = np.linalg.norm(np.dot(self.nodes[self.TSnode].gradient.T,self.nodes[self.TSnode].constraints[:,0])*self.nodes[self.TSnode].constraints[:,0])
print(" ts_cgradq %5.4f" % ts_cgradq)
ts_gradrms=self.nodes[self.TSnode].gradrms
self.dE_iter=abs(self.emax-self.emaxp)
print(" dE_iter ={:2.2f}".format(self.dE_iter))
# => calculate totalgrad <= #
totalgrad,gradrms,sum_gradrms = self.calc_optimization_metrics(self.nodes)
# => Check Convergence <= #
if self.is_converged(totalgrad,fp,rtype,ts_cgradq):
break
# => set stage <= #
form_TS_hess = self.set_stage(totalgrad,sum_gradrms,ts_cgradq,ts_gradrms,fp)
# => Reparam the String <= #
if oi!=max_iter-1:
self.reparameterize(nconstraints=nconstraints)
# store reparam energies
print(" V_profile (after reparam): ", end=' ')
self.print_energies()
## Modify TS Hess if necessary ##
# from set stage
if form_TS_hess:
self.ictan,self.dqmaga = self.get_three_way_tangents(self.nodes,self.energies)
self.modify_TS_Hess()
if self.optimizer[self.TSnode].options['DMAX']>0.1:
self.optimizer[self.TSnode].options['DMAX']=0.1
# opt decided Hess is not good because of overlap
elif self.find and not self.optimizer[self.TSnode].maxol_good:
self.ictan,self.dqmaga = self.get_three_way_tangents(self.nodes,self.energies)
self.modify_TS_Hess()
elif self.find and (self.optimizer[self.TSnode].nneg > 3 or self.optimizer[self.TSnode].nneg==0 or self.hess_counter > 3 or (self.TS_E_0 - self.emax) > 10.) and ts_gradrms >self.CONV_TOL:
if self.hessrcount<1 and self.pTSnode == self.TSnode:
print(" resetting TS node coords Ut (and Hessian)")
self.ictan,self.dqmaga = self.get_three_way_tangents(self.nodes,self.energies)
self.modify_TS_Hess()
self.nhessreset=10
self.hessrcount=1
else:
print(" Hessian consistently bad, going back to climb (for 3 iterations)")
self.find=False
#self.optimizer[self.TSnode] = beales_cg(self.optimizer[self.TSnode].options.copy().set_values({"Linesearch":"backtrack"}))
self.nclimb=2
#elif self.find and self.optimizer[self.TSnode].nneg > 1 and ts_gradrms < self.CONV_TOL:
# print(" nneg > 1 and close to converging -- reforming Hessian")
# self.ictan,self.dqmaga = self.get_three_way_tangents(self.nodes,self.energies)
# self.modify_TS_Hess(self.TSnode)
elif self.find and self.optimizer[self.TSnode].nneg <= 3:
self.hessrcount-=1
self.hess_counter += 1
if self.pTSnode!=self.TSnode and self.climb:
#self.optimizer[self.TSnode] = beales_cg(self.optimizer[self.TSnode].options.copy())
#self.optimizer[self.pTSnode] = self.optimizer[0].__class__(self.optimizer[self.TSnode].options.copy())
if self.climb and not self.find:
print(" slowing down climb optimization")
self.optimizer[self.TSnode].options['DMAX'] /= self.newclimbscale
self.optimizer[self.TSnode].options['SCALEQN'] = 2.
if self.optimizer[self.TSnode].SCALE_CLIMB <5.:
self.optimizer[self.TSnode].SCALE_CLIMB +=1.
self.optimizer[self.pTSnode].options['SCALEQN'] = 1.
self.ts_exsteps=1
if self.newclimbscale<5.0:
self.newclimbscale +=1.
elif self.find:
self.find = False
self.nclimb=1
print(" Find bad, going back to climb")
#self.optimizer[self.TSnode] = beales_cg(self.optimizer[self.pTSnode].options.copy().set_values({"Linesearch":"backtrack"}))
#self.optimizer[self.pTSnode] = self.optimizer[0].__class__(self.optimizer[self.TSnode].options.copy())
#self.self.ictan,self.dqmaga = self.get_three_way_tangents(self.nodes,self.energies)
#self.modify_TS_Hess(self.TSnode)
# => write Convergence to file <= #
filename = 'scratch/opt_iters_{:03}_{:03}.xyz'.format(self.ID,oi)
write_molden_geoms(filename,self.geometries,self.energies,self.gradrmss,self.dEs)
#TODO prints tgrads and jobGradCount
print("opt_iter: {:2} totalgrad: {:4.3} gradrms: {:5.4} max E({}) {:5.4}".format(oi,float(totalgrad),float(gradrms),self.TSnode,float(self.emax)))
print('\n')
#TODO Optimize TS node to a finer convergence
#if rtype==2:
filename="opt_converged_{:03d}.xyz".format(self.ID)
print(" Printing string to " + filename)
write_molden_geoms(filename,self.geometries,self.energies,self.gradrms,self.dEs)
return
print_energies(self)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def print_energies(self):
for n in range(len(self.energies)):
print(" {:7.3f}".format(float(self.energies[n])), end=' ')
print()
refresh_coordinates(self, update_TS = False)
inherited
🔗
Refresh the DLC coordinates for the string
Source code in pygsm/growing_string_methods/de_gsm.py
def refresh_coordinates(self,update_TS=False):
'''
Refresh the DLC coordinates for the string
'''
energies = self.energies
TSnode = self.TSnode
if not self.done_growing:
#TODO
for n in range(1,self.nnodes-1):
if self.nodes[n] is not None:
self.nodes[n].update_coordinate_basis(self.ictan[n])
else:
for n in range(1,self.nnodes-1):
# don't update tsnode coord basis
if n!=TSnode or (n==TSnode and update_TS):
self.nodes[n].update_coordinate_basis(self.ictan[n])
# update newic coordinate basis
#self.newic.xyz = self.nodes[n].xyz
#Vecs = self.newic.update_coordinate_basis(self.ictan[n])
#self.nodes[n].coord_basis = Vecs
reparameterize(self, ic_reparam_steps = 8, n0 = 0, nconstraints = 1, rtype = 0)
inherited
🔗
Reparameterize the string
Source code in pygsm/growing_string_methods/de_gsm.py
def reparameterize(self,ic_reparam_steps=8,n0=0,nconstraints=1,rtype=0):
'''
Reparameterize the string
'''
print(self.interp_method)
if self.interp_method == 'DLC':
print('reparameterizing')
self.ic_reparam(nconstraints=nconstraints)
elif self.interp_method == 'Geodesic':
self.geodesic_reparam()
return
set_active(self, nR, nP)
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def set_active(self,nR,nP):
#print(" Here is active:",self.active)
if nR!=nP and self.growth_direction==0:
print((" setting active nodes to %i and %i"%(nR,nP)))
elif self.growth_direction==1:
print((" setting active node to %i "%nR))
elif self.growth_direction==2:
print((" setting active node to %i "%nP))
else:
print((" setting active node to %i "%nR))
for i in range(self.nnodes):
if self.nodes[i] != None:
self.optimizer[i].conv_grms = self.options['CONV_TOL']*2.
self.optimizer[nR].conv_grms = self.options['ADD_NODE_TOL']
self.optimizer[nP].conv_grms = self.options['ADD_NODE_TOL']
print(" conv_tol of node %d is %.4f" % (nR,self.optimizer[nR].conv_grms))
print(" conv_tol of node %d is %.4f" % (nP,self.optimizer[nP].conv_grms))
self.active[nR] = True
self.active[nP] = True
if self.growth_direction==1:
self.active[nP]=False
if self.growth_direction==2:
self.active[nR]=False
#print(" Here is new active:",self.active)
set_finder(self, rtype)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def set_finder(self,rtype):
assert rtype in [0,1,2], "rtype not defined"
print('')
print("*********************************************************************")
if rtype==2:
print("****************** set climber and finder to True *******************")
self.climber=True
self.finder=True
elif rtype==1:
print("***************** setting climber to True*************************")
self.climber=True
else:
print("******** Turning off climbing image and exact TS search **********")
print("*********************************************************************")
set_opt_type(self, n, quiet = False)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def set_opt_type(self,n,quiet=False):
#TODO error for seam climb
opt_type='ICTAN'
if self.climb and n==self.TSnode and not self.find and self.nodes[n].PES.__class__.__name__!="Avg_PES":
opt_type='CLIMB'
#opt_type='BEALES_CG'
elif self.find and n==self.TSnode:
opt_type='TS'
elif self.nodes[n].PES.__class__.__name__=="Avg_PES":
opt_type='SEAM'
if self.climb and n==self.TSnode:
opt_type='TS-SEAM'
if not quiet:
print((" setting node %i opt_type to %s" %(n,opt_type)))
#if isinstance(self.optimizer[n],beales_cg) and opt_type!="BEALES_CG":
# raise RuntimeError("This shouldn't happen")
return opt_type
set_stage(self, totalgrad, sumgradrms, ts_cgradq, ts_gradrms, fp)
inherited
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def set_stage(self,totalgrad,sumgradrms, ts_cgradq,ts_gradrms,fp):
form_TS_hess=False
sum_conv_tol = np.sum([self.optimizer[n].conv_grms for n in range(1,self.nnodes-1)])+ 0.0005
#TODO totalgrad is not a good criteria for large systems
if (totalgrad < 0.15 or sumgradrms<sum_conv_tol or ts_cgradq < 0.01) and fp>0 and self.dE_iter < 1.: # extra criterion in og-gsm for added
if not self.climb and self.climber:
print(" ** starting climb **")
self.climb=True
print(" totalgrad %5.4f gradrms: %5.4f gts: %5.4f" %(totalgrad,ts_gradrms,ts_cgradq))
# set to beales
#self.optimizer[self.TSnode] = beales_cg(self.optimizer[self.TSnode].options.copy().set_values({"Linesearch":"backtrack"}))
#self.optimizer[self.TSnode].options['DMAX'] /= self.newclimbscale
# overwrite this here just in case TSnode changed wont cause slow down climb
self.pTSnode = self.TSnode
elif (self.climb and not self.find and self.finder and self.nclimb<1 and
((totalgrad<0.2 and ts_gradrms<self.options['CONV_TOL']*10. and ts_cgradq<0.01) or #
(totalgrad<0.1 and ts_gradrms<self.options['CONV_TOL']*10. and ts_cgradq<0.02) or #
(sumgradrms< sum_conv_tol) or
(ts_gradrms<self.options['CONV_TOL']*5.) # used to be 5
)) and self.dE_iter<1. :
print(" ** starting exact climb **")
print(" totalgrad %5.4f gradrms: %5.4f gts: %5.4f" %(totalgrad,ts_gradrms,ts_cgradq))
self.find=True
form_TS_hess=True
self.optimizer[self.TSnode] = eigenvector_follow(self.optimizer[self.TSnode].options.copy())
print(type(self.optimizer[self.TSnode]))
self.optimizer[self.TSnode].options['SCALEQN'] = 1.
self.nhessreset=10 # are these used??? TODO
self.hessrcount=0 # are these used?! TODO
if self.climb:
self.nclimb-=1
#for n in range(1,self.nnodes-1):
# self.active[n]=True
# self.optimizer[n].options['OPTTHRESH']=self.options['CONV_TOL']*2
self.nhessreset-=1
return form_TS_hess
set_V0(self)
🔗
Source code in pygsm/growing_string_methods/de_gsm.py
def set_V0(self):
self.nodes[0].V0 = self.nodes[0].energy
#TODO should be actual gradient
self.nodes[0].gradrms = 0.
if self.growth_direction!=1:
self.nodes[-1].gradrms = 0.
print(" Energy of the end points are %4.3f, %4.3f" %(self.nodes[0].energy,self.nodes[-1].energy))
print(" relative E %4.3f, %4.3f" %(0.0,self.nodes[-1].energy-self.nodes[0].energy))
else:
print(" Energy of end points are %4.3f " % self.nodes[0].energy)
#self.nodes[-1].energy = self.nodes[0].energy
#self.nodes[-1].gradrms = 0.
setup_from_geometries(self, input_geoms, reparametrize = False, restart_energies = True)
inherited
🔗
Restart
Source code in pygsm/growing_string_methods/de_gsm.py
def setup_from_geometries(self,input_geoms,reparametrize=False,restart_energies=True):
'''
Restart
'''
printcool("Restarting GSM from geometries")
self.growth_direction=0
nstructs=len(input_geoms)
if nstructs != self.nnodes:
print('need to interpolate')
#if self.interp_method=="DLC": TODO
symbols = get_atoms(input_geoms[0])
old_xyzs = [ xyz_to_np( geom ) for geom in input_geoms ]
xyzs = redistribute(symbols,old_xyzs,self.nnodes,tol=2e-3*5)
geoms = [ np_to_xyz(input_geoms[0],xyz) for xyz in xyzs ]
nstructs = len(geoms)
else:
geoms = input_geoms
self.gradrms = [0.]*nstructs
self.dE = [1000.]*nstructs
self.isRestarted=True
self.done_growing=True
# set coordinates from geoms
self.nodes[0].xyz = xyz_to_np(geoms[0])
self.nodes[nstructs-1].xyz = xyz_to_np(geoms[-1])
for struct in range(1,nstructs-1):
self.nodes[struct] = Molecule.copy_from_options(self.nodes[struct-1],
xyz_to_np(geoms[struct]),
new_node_id=struct,
copy_wavefunction=False)
self.nodes[struct].newHess=5
# Turning this off
#self.nodes[struct].gradrms = np.sqrt(np.dot(self.nodes[struct].gradient,self.nodes
#self.nodes[struct].gradrms=grmss[struct]
#self.nodes[struct].PES.dE = dE[struct]
self.nnodes=self.nR=nstructs
if reparametrize:
printcool("Reparametrizing")
self.reparameterize(ic_reparam_steps=8)
write_xyz_files('grown_string1_{:03}.xyz'.format(self.ID))
if restart_energies:
# initial energy
self.nodes[0].V0 = self.nodes[0].energy
self.energies[0] = 0.
print(" initial energy is %3.4f" % self.nodes[0].energy)
for struct in range(1,nstructs-1):
print(" energy of node %i is %5.4f" % (struct,self.nodes[struct].energy))
self.energies[struct] = self.nodes[struct].energy - self.nodes[0].V0
print(" Relative energy of node %i is %5.4f" % (struct,self.energies[struct]))
print(" V_profile: ", end=' ')
energies= self.energies
for n in range(self.nnodes):
print(" {:7.3f}".format(float(energies[n])), end=' ')
print()
print(" setting all interior nodes to active")
for n in range(1,self.nnodes-1):
self.active[n]=True
self.optimizer[n].conv_grms=self.options['CONV_TOL']*2.5
self.optimizer[n].options['DMAX'] = 0.05
return