- example script for explicit solvation with counterion

- calc_hessian_xtb: choose between xtbmethods
- define_topology and write_pdbfile_openmm (resname can be chosen)
- insert_solute_into_solvent: now 2 solutes can be inserted at the same time, convenient for counterions
- geometric: re-enabled first and each Hessian options for TSOpt jobs
- NEBTS: Enabled single-iter NEBTS with partial Hessian. Added keywords: partial_hessian_atoms and partial_npoint_choice

Author: RagnarB83

ash/interfaces/interface_geometric_new.py

@@ -383,13 +383,15 @@ def hessian_option(self,fragment,actatoms,theory,charge,mult,modelhessian):
                 elif "file:" in self.hessian:
                     hessianfile = self.hessian.replace("file:","")
 
-                print("Checking that defined Hessian is compatible with active-region")
-                hessian_read = read_hessian(hessianfile)
-                print("actatoms:", actatoms)
-                if hessian_read.shape[0] != 3*len(atomsused):
-                    print(f"Error: Hessian shape is {hessian_read.shape}  which is incompatible with the  number of active atoms present ({len(atomsused)})")
-                    print(f"Hessian should have dimension of 3*N x 3*N where N is the number of active-atoms of the system (should be : {3*len(atomsused)} x {3*len(atomsused)})")
-                    ashexit()
+                #Allow first and each options still
+                if self.hessian not in ['first','each']:
+                    print("Checking that defined Hessian is compatible with active-region")
+                    hessian_read = read_hessian(hessianfile)
+                    print("actatoms:", actatoms)
+                    if hessian_read.shape[0] != 3*len(atomsused):
+                        print(f"Error: Hessian shape is {hessian_read.shape}  which is incompatible with the  number of active atoms present ({len(atomsused)})")
+                        print(f"Hessian should have dimension of 3*N x 3*N where N is the number of active-atoms of the system (should be : {3*len(atomsused)} x {3*len(atomsused)})")
+                        ashexit()
 
             elif self.hessian == None:
                 if self.printlevel >= 1:

ash/interfaces/interface_knarr.py

@@ -77,7 +77,8 @@ def NEBTS(reactant=None, product=None, theory=None, images=8, CI=True, free_end=
         tol_turn_on_ci=1.0,  runmode='serial', numcores=1, charge=None, mult=None, printlevel=1, ActiveRegion=False, actatoms=None,
         interpolation="IDPP", idpp_maxiter=700, idpp_springconst=5.0, restart_file=None, TS_guess_file=None, mofilesdir=None,
         OptTS_maxiter=100, OptTS_print_atoms_list=None, OptTS_convergence_setting=None, OptTS_conv_criteria=None, OptTS_coordsystem='tric',
-        hessian_for_TS=None, modelhessian='unit', tsmode_tangent_threshold=0.1, subfrctor=1):
+        hessian_for_TS=None, modelhessian='unit', tsmode_tangent_threshold=0.1, subfrctor=1,
+        partial_hessian_atoms=None, partial_npoint_choice=2):
 
     print_line_with_mainheader("NEB+TS")
     module_init_time=time.time()
@@ -181,8 +182,22 @@ def NEBTS(reactant=None, product=None, theory=None, images=8, CI=True, free_end=
     elif hessian_for_TS == 'each':
         hessianoption="each"
     #xTB Hessian option
-    elif hessian_for_TS == 'xtb':
-        hessianfile = calc_hessian_xtb(fragment=SP, runmode='serial', actatoms=actatoms, numcores=max_cores, use_xtb_feature=True)
+    elif 'xtb' in hessian_for_TS.lower():
+        print("xtB-type Hessian chosen (xtb found in string)")
+        if 'gfn1' in hessian_for_TS.lower():
+            print("Using GFN1-xTB Hessian")
+            xtbmethod="GFN1"
+        elif 'gfn2' in hessian_for_TS.lower():
+            print("Using GFN2-xTB Hessian")
+            xtbmethod="GFN2"
+        else:
+            print("Neither gfn1 or gfn2 was chosen (control this by choosing hessian_for_TS='GFN1-xTB' or hessian_for_TS='GFN2-xTB' )")
+            print("Choosing xtbmethod to be GFN2-xTB by default")
+            xtbmethod="GFN2"
+
+
+        hessianfile = calc_hessian_xtb(fragment=SP, runmode='serial', actatoms=actatoms, numcores=max_cores, use_xtb_feature=True,
+                                       xtbmethod=xtbmethod)
         hessianoption='file:'+str(hessianfile)
     #Cheap model Hessian
     #NOTE: None of these work well.  Need to use tangent to modify
@@ -207,35 +222,54 @@ def NEBTS(reactant=None, product=None, theory=None, images=8, CI=True, free_end=
     # Add connecting atoms and erform partial Hessian optimization
     elif hessian_for_TS == 'partial':
         print("hessian_for_TS option: partial")
-        print(f"Using climbing image tangent to find dominant atoms in approximate TS mode (tsmode_tangent_threshold={tsmode_tangent_threshold})")
-
-        #Getting tangent and atoms that contribute at least X to tangent where X is tsmode_tangent_threshold=0.1 (default)
-        tangent = read_tangent("CItangent.xyz")
-        TSmodeatoms = list(np.where(np.any(abs(tangent)>tsmode_tangent_threshold, axis=1))[0])
-        #Convert activeregion atom indices to full system indices
-        print("Determining atoms contributing the most to TS mode")
-        if ActiveRegion is True:
-            print("TSmodeatoms (active region):", TSmodeatoms)
-            TSmodeatoms = [actatoms[a] for a in TSmodeatoms]
-            print("TSmodeatoms (full system):", TSmodeatoms)
+        print("Using partial Hessian to approximate Hessian for saddle-point optimization.")
+
+
+        if Singleiter is True:
+            print("This is a Singleiter NEBTS job.")
+            if partial_hessian_atoms is None:
+                print("Error: partial_hessian_atoms is not set")
+                print("In Singleiter NEBTS jobs no climbing image tangent is available. Setting the partial_hessian_atoms list is necessary.")
+                print("Example: partial_hessian_atoms=[15,16,18,19]  where the numbers are atom indices believed to contribute the most to the TS mode")
+                print("Exiting")
+                ashexit()
+            Final_partial_hessatoms=partial_hessian_atoms
+            print(f"Performing partial Hessian calculation using atom-list: {Final_partial_hessatoms}")
         else:
-            print("TSmodeatoms (full system):", TSmodeatoms)
-        #Now finding the atoms that TSmodeatoms are connected to, for both R, P and SP
-        print("Now finding connected atoms to TSmode-atoms")
-        result_R = get_conn_atoms_for_list(fragment=reactant, atoms=TSmodeatoms)
-        print("result_R:", result_R)
-        result_P = get_conn_atoms_for_list(fragment=product, atoms=TSmodeatoms)
-        print("result_P:", result_P)
-        result_SP = get_conn_atoms_for_list(fragment=SP, atoms=TSmodeatoms)
-        print("result_SP:", result_SP)
-        #Combining
-        Final_partial_hessatoms = np.unique(result_R + result_P + result_SP).tolist()
-
-        print(f"Performing partial Hessian calculation using atom-list: {Final_partial_hessatoms}")
-        print("This corresponds to atoms contributing to the TS-mode and connected atoms")
+            print("This is a regular NEBTS job. Partial Hessian can either be approximated from tangent or explicity set via partial_hessian_atoms keyword")
+            if partial_hessian_atoms is None:
+                print("partial_hessian_atoms list is empty. Will now try to use climbing image tangent to guess the important atoms")
+                print(f"Using climbing image tangent to find dominant atoms in approximate TS mode (tsmode_tangent_threshold={tsmode_tangent_threshold})")
+
+                #Getting tangent and atoms that contribute at least X to tangent where X is tsmode_tangent_threshold=0.1 (default)
+                tangent = read_tangent("CItangent.xyz")
+                TSmodeatoms = list(np.where(np.any(abs(tangent)>tsmode_tangent_threshold, axis=1))[0])
+                #Convert activeregion atom indices to full system indices
+                print("Determining atoms contributing the most to TS mode")
+                if ActiveRegion is True:
+                    print("TSmodeatoms (active region):", TSmodeatoms)
+                    TSmodeatoms = [actatoms[a] for a in TSmodeatoms]
+                    print("TSmodeatoms (full system):", TSmodeatoms)
+                else:
+                    print("TSmodeatoms (full system):", TSmodeatoms)
+                #Now finding the atoms that TSmodeatoms are connected to, for both R, P and SP
+                print("Now finding connected atoms to TSmode-atoms")
+                result_R = get_conn_atoms_for_list(fragment=reactant, atoms=TSmodeatoms)
+                print("result_R:", result_R)
+                result_P = get_conn_atoms_for_list(fragment=product, atoms=TSmodeatoms)
+                print("result_P:", result_P)
+                result_SP = get_conn_atoms_for_list(fragment=SP, atoms=TSmodeatoms)
+                print("result_SP:", result_SP)
+                #Combining
+                Final_partial_hessatoms = np.unique(result_R + result_P + result_SP).tolist()
+                print(f"Performing partial Hessian calculation using atom-list: {Final_partial_hessatoms}")
+                print("This corresponds to atoms contributing to the TS-mode and connected atoms")
+            else:
+                Final_partial_hessatoms=partial_hessian_atoms
+                print(f"Performing partial Hessian calculation using atom-list: {Final_partial_hessatoms}")
         #TODO: Option to run this in parallel ?
         #Or just enable theory parallelization
-        result_freq = ash.NumFreq(theory=theory, fragment=SP, printlevel=0, npoint=2, hessatoms=Final_partial_hessatoms, runmode=runmode, numcores=numcores)
+        result_freq = ash.NumFreq(theory=theory, fragment=SP, printlevel=0, npoint=partial_npoint_choice, hessatoms=Final_partial_hessatoms, runmode=runmode, numcores=numcores)
 
         #Combine partial exact Hessian with model Hessian(Almloef, Lindh, Schlegel or unit)
         #Large Hessian is the actatoms Hessian if actatoms provided

ash/modules/module_coords.py

@@ -760,7 +760,7 @@ def write_pdbfile(self,filename="Fragment"):
         return f"{filename}.pdb"
 
     # Create new topology from scratch if none is defined (defined automatically when reading PDB-files by OpenMM)
-    def define_topology(self, scale=1.0, tol=0.1):
+    def define_topology(self, scale=1.0, tol=0.1, resname="MOL"):
         try:
             import openmm.app
         except ImportError:
@@ -769,7 +769,7 @@ def define_topology(self, scale=1.0, tol=0.1):
         print("Defining new basic single-chain, single-residue topology")
         self.pdb_topology = openmm.app.Topology()
         chain = self.pdb_topology.addChain()
-        residue = self.pdb_topology.addResidue("MOL", chain)
+        residue = self.pdb_topology.addResidue(resname, chain)
 
         # Defaultdictionary to keep track of unique element-atomnames
         atomnames_dict=defaultdict(int)
@@ -788,7 +788,7 @@ def define_topology(self, scale=1.0, tol=0.1):
 
     # Write PDB-file via OpenMM
     def write_pdbfile_openmm(self,filename="Fragment", calc_connectivity=False, pdb_topology=None,
-                             skip_connectivity=False):
+                             skip_connectivity=False, resname="MOL"):
         print("write_pdbfile_openmm\n")
         try:
             import openmm.app
@@ -807,7 +807,7 @@ def write_pdbfile_openmm(self,filename="Fragment", calc_connectivity=False, pdb_
             print("Warning: ASH Fragment has no PDB-file topology defined (required for PDB-file writing)")
             print("Now defining new topology from scratch")
             if pdb_topology is None:
-                self.define_topology() #Creates self.pdb_topology
+                self.define_topology(resname=resname) #Creates self.pdb_topology
         else:
             print("Using pdbtopology found in ASH fragment")
 
@@ -4043,13 +4043,29 @@ def bounding_box_dimensions(coordinates,shift=0.0):
     final_dims = dimensions + shift
     return dimensions  # Return the dimensions of the bounding box
 
+# Combien and place 2 fragments
+def combine_and_place_fragments(ref_frag, trans_frag):
+    
+    for displacement in [0.0, 2.0, 4.0, 6.0,8.0,10.0,12.0,14.0]:
+        #Translating 2nd frag by displacement in 
+        trans_frag.coords[:,-1] +=displacement
+        members = get_molecule_members_loop_np2(np.vstack((ref_frag.coords,trans_frag.coords)), ref_frag.elems+trans_frag.elems, 
+                                        10, 1.0, 0.4, atomindex=0, membs=None)
+        if len(members) == ref_frag.numatoms:
+            print("Molecules are sufficiently far apart")
+            break
+    
+    combined_solute = Fragment(elems=ref_frag.elems+trans_frag.elems, coords = np.vstack((ref_frag.coords,trans_frag.coords)), 
+                               printlevel=2)
+    return combined_solute
+
 
 #Simple function to combine 2 ASH fragments where one is assumed to be a solute (fewer atoms) and the other assumed to be
 #some kind of solvent system (box,sphere etc.)
 #Use tolerance (tol) e.g. to control how many solvent molecules around get deleted
 #Currently using 0.4 as default based on threonine in acetonitrile example
-def insert_solute_into_solvent(solute=None, solvent=None, scale=1.0, tol=0.4, write_pdb=False, write_solute_connectivity=True,
-                                       solute_pdb=None, solvent_pdb=None, outputname="solution.pdb", write_PBC_info=True):
+def insert_solute_into_solvent(solute=None, solute2=None, solvent=None, scale=1.0, tol=0.4, write_pdb=False, write_solute_connectivity=True,
+                                       solute_pdb=None, solute2_pdb=None, solvent_pdb=None, outputname="solution.pdb", write_PBC_info=True):
     print("\ninsert_solute_into_solvent\n")
     #Early exits
     if write_pdb:
@@ -4060,32 +4076,63 @@ def insert_solute_into_solvent(solute=None, solvent=None, scale=1.0, tol=0.4, wr
     if solute is None and solute_pdb is not None:
         print("No solute fragment provided but solute_pdb is set. Reading solute fragment from PDB-file")
         solute = Fragment(pdbfile=solute_pdb)
+    if solute2 is None and solute2_pdb is not None:
+        print("No solute2 fragment provided but solute2_pdb is set. Reading solute2 fragment from PDB-file")
+        solute2 = Fragment(pdbfile=solute2_pdb)
     if solvent is None and solvent_pdb is not None:
         print("No solvent fragment provided but solvent_pdb is set. Reading solvent fragment from PDB-file")
         solvent = Fragment(pdbfile=solvent_pdb)
 
     #Get centers
     com_box = solvent.get_coordinate_center()
-    com_solute = solute.get_coordinate_center()
 
-    #Translating solute coords
-    trans_coord=np.array(com_box)-np.array(com_solute)
-    solute.coords=solute.coords+trans_coord
+    if solute2 is  None:
+        com_solute = solute.get_coordinate_center()
+        #Translating solute coords
+        trans_coord=np.array(com_box)-np.array(com_solute)
+        solute.coords=solute.coords+trans_coord
+    else:
+        #Combine and Translate solute2 so that it is close to solute1
+        combined_solute = combine_and_place_fragments(ref_frag=solute, trans_frag=solute2)
+
+        #COM and trans
+        com_solute = combined_solute.get_coordinate_center()
+        trans_coord=np.array(com_box)-np.array(com_solute)
+        #Translate combined solute fragment
+        combined_solute.coords=combined_solute.coords+trans_coord
+
 
     #New Fragment by combining element lists and coordinates
-    new_frag = Fragment(elems=solute.elems+solvent.elems, coords = np.vstack((solute.coords,solvent.coords)), printlevel=2)
+    if solute2 is None:
+        print("Combining solute and solvent")
+        new_frag = Fragment(elems=solute.elems+solvent.elems, coords = np.vstack((solute.coords,solvent.coords)), printlevel=2)
+    else:
+        print("Combining combined_solute and solvent")
+        new_frag = Fragment(elems=combined_solute.elems+solvent.elems, coords = np.vstack((combined_solute.coords,solvent.coords)), printlevel=2)
 
     new_frag.write_xyzfile(xyzfilename="solution-pre.xyz")
     #Trim by removing clashing atoms
     new_frag.printlevel=1
+
     #Find atoms connected to solute index 0. Uses scale and tol
     membs = get_molecule_members_loop_np2(new_frag.coords, new_frag.elems, 20, scale, tol, atomindex=0, membs=None)
-    print("Found clashing solvent atoms:", membs)
     delatoms=[]
     for i in membs:
         if i >= solute.numatoms:
             delatoms.append(i)
+    print("First delatoms:", delatoms)
+    if solute2 is not None:
+        membs2 = get_molecule_members_loop_np2(new_frag.coords, new_frag.elems, 20, scale, tol, atomindex=solute.numatoms, membs=None)
+        print("membs2:", membs2)
+        for j in membs2:
+            if j >= solute.numatoms+solute2.numatoms:
+                if j not in delatoms:
+                    delatoms.append(j)
+    print("Final delatoms:", delatoms)
+    
+    # Deleting
     delatoms.sort(reverse=True)
+    print("Found clashing solvent atoms:", delatoms)
     for d in delatoms:
         new_frag.delete_atom(d)
     new_frag.printlevel=2
@@ -4113,7 +4160,18 @@ def insert_solute_into_solvent(solute=None, solvent=None, scale=1.0, tol=0.4, wr
 
         #Create modeller object
         modeller = openmm.app.Modeller(pdb1.topology, pdb1.positions) #Add pdbfile1
+        
+        #solute2
+        if solute2 is not None:
+            print("Adding solute2")
+            pdb_solute2 = openmm.app.PDBFile(solute2_pdb)
+            solute2_resname= list(pdb_solute2.topology.residues())[0].name
+            print("solute2_resname:", solute2_resname)
+            modeller.add(pdb_solute2.topology, pdb_solute2.positions) #Add pdbfile2
+        print("Adding solvent")
         modeller.add(pdb2.topology, pdb2.positions) #Add pdbfile2
+
+
         #Delete clashing atoms from topology
         toDelete = [r for j,r in enumerate(modeller.topology.atoms()) if j in delatoms]
         modeller.delete(toDelete)

ash/modules/module_freq.py

@@ -1874,7 +1874,7 @@ def read_hessian(file):
 
 #Calculate Hessian (GFN1) for a fragment.
 def calc_hessian_xtb(fragment=None, runmode='serial', actatoms=None, numcores=1, use_xtb_feature=True,
-    charge=None, mult=None):
+    charge=None, mult=None, xtbmethod="GFN1"):
 
     print_line_with_mainheader("calc_hessian_xtb")
 
@@ -1890,7 +1890,7 @@ def calc_hessian_xtb(fragment=None, runmode='serial', actatoms=None, numcores=1,
         fragment = ash.Fragment(elems=subelems,coords=subcoords, printlevel=0)
     print("Will now calculate xTB Hessian")
     #Creating xtb theory object
-    xtb = ash.xTBTheory(xtbmethod='GFN1', numcores=numcores)
+    xtb = ash.xTBTheory(xtbmethod=xtbmethod, numcores=numcores)
 
     #Get Hessian from xTB directly (avoiding ASH NumFreq)
     if use_xtb_feature == True: