Read bond types from CIF file.

.vscode/settings.json

@@ -0,0 +1,4 @@
+{
+    "python-envs.defaultEnvManager": "ms-python.python:conda",
+    "python-envs.defaultPackageManager": "ms-python.python:conda"
+}

lammps_interface/ForceFields.py

@@ -2150,6 +2150,8 @@ def bond_term(self, edge):
         if (self.keep_metal_geometry) and (n1_data['atomic_number'] in METALS
             or n2_data['atomic_number'] in METALS):
             r0 = data['length']
+            K *= 50
+
         data['potential'] = BondPotential.Harmonic()
         data['potential'].K = K
         data['potential'].R0 = r0
@@ -2228,6 +2230,19 @@ def angle_term(self, angle):
         ka = beta*(za*zc /(r_ac**5.))*r_ab*r_bc
         ka *= (3.*r_ab*r_bc*(1. - cosT0*cosT0) - r_ac*r_ac*cosT0)
 
+        if self.keep_metal_geometry:
+            is_a_metal = a_data['atomic_number'] in METALS
+            is_b_metal = b_data['atomic_number'] in METALS
+            is_c_metal = c_data['atomic_number'] in METALS
+
+            # If ANY atom in the angle is a metal, lock it with a massive harmonic constant
+            if any([is_a_metal, is_b_metal, is_c_metal]):
+                theta0 = self.graph.compute_angle_between(a, b, c)
+                data['potential'] = AnglePotential.Harmonic()
+                data['potential'].K = max(500.0, ka)  # Masssive stiffness acts as a rigid constraint
+                data['potential'].theta0 = theta0
+                return 1
+
         if angle_type in sf or (angle_type == 'tetrahedral' and int(theta0) == 90):
             if angle_type == 'linear':
                 kappa = ka
@@ -2499,8 +2514,12 @@ def detect_ff_terms(self):
                             if neigh_elem <= set(["Si", "Al"]):
                                 data['force_field_type'] = "O_3_z"
                         elif data['element'] == "S":
-                            # default sp3 hybridized sulphur set to S_3+6
-                            data['force_field_type'] = "S_3+6"
+                            if self.graph.degree(node) == 4:
+                                data['force_field_type'] = "S_3+6"
+                            elif self.graph.degree(node) == 3:
+                                data['force_field_type'] = "S_3+4"
+                            elif self.graph.degree(node) == 2:
+                                data['force_field_type'] = "S_3+2"
 
                     elif data['hybridization'] == "aromatic":
                         data['force_field_type'] = "%s_R"%data['element']
@@ -3175,6 +3194,8 @@ def bond_term(self, edge):
         if (self.keep_metal_geometry) and (n1_data['atomic_number'] in METALS
             or n2_data['atomic_number'] in METALS):
             r0 = data['length']
+            K *= 50
+
         data['potential'] = BondPotential.Harmonic()
         data['potential'].K = K
         data['potential'].R0 = r0
@@ -3247,6 +3268,20 @@ def angle_term(self, angle):
         beta = 664.12/r_ab/r_bc
         ka = beta*(za*zc /(r_ac**5.))*r_ab*r_bc
         ka *= (3.*r_ab*r_bc*(1. - cosT0*cosT0) - r_ac*r_ac*cosT0)
+
+        if self.keep_metal_geometry:
+            is_a_metal = a_data['atomic_number'] in METALS
+            is_b_metal = b_data['atomic_number'] in METALS
+            is_c_metal = c_data['atomic_number'] in METALS
+
+            # If ANY atom in the angle is a metal, lock it with a massive harmonic constant
+            if any([is_a_metal, is_b_metal, is_c_metal]):
+                theta0 = self.graph.compute_angle_between(a, b, c)
+                data['potential'] = AnglePotential.Harmonic()
+                data['potential'].K = max(500.0, ka)  # Masssive stiffness acts as a rigid constraint
+                data['potential'].theta0 = theta0
+                return 1
+
         #if ("special_flag" in b_data.keys()) and b_data["special_flag"] == "Cu_pdw":
         #    angle_type = "None"
         #    print(self.graph.compute_angle_between(a, b, c))
@@ -3297,9 +3332,10 @@ def angle_term(self, angle):
             kappa = ka
             data['potential'] = AnglePotential.Fourier()
             data['potential'].K = kappa
-            data['potential'].C0 = c0
-            data['potential'].C1 = c1
-            data['potential'].C2 = c2
+            data['potential'].C0 = c0 if abs(c0) < 100 else 0
+            data['potential'].C1 = c1 if abs(c1) < 100 else 0
+            data['potential'].C2 = c2 if abs(c2) < 100 else 0
+
             return 1
 
     def uff_angle_type(self, b):
@@ -3616,6 +3652,14 @@ def detect_ff_terms(self):
                         data['force_field_type'] = "%s_2"%data['element']
                     elif data['hybridization'] == "sp":
                         data['force_field_type'] = "%s_1"%data['element']
+
+                    if data['element'] == "S":
+                        if self.graph.degree(node) == 4:
+                            data['force_field_type'] = "S_3+6"
+                        elif self.graph.degree(node) == 3:
+                            data['force_field_type'] = "S_3+4"
+                        elif self.graph.degree(node) == 2:
+                            data['force_field_type'] = "S_3+2"
                 if data['element'] == "O" and self.graph.degree(node) == 2:
                     if neigh_elem <= metals:
                         data['force_field_type'] = "O_2"

lammps_interface/InputHandler.py

@@ -239,6 +239,14 @@ def run_command_line_options(self):
                                      help="Max deviation of adjusted variable "+
                                           "at each step is scaled by MAX_DEV/ITER_COUNT. "+
                                           "Default is 0.01 (ideal for volume).")
+        parameter_group.add_argument("--min_eval",
+                                     action="store",
+                                     type=float,
+                                     default=0.000001,
+                                     dest="min_eval",
+                                     help="Minimal change in energy per atom to consider" +
+                                          "the optimization converged." +
+                                          "Default is 1e-6 (ideal for most cases).")
         parameter_group.add_argument("--temperature",
                                      action="store",
                                      type=float,

lammps_interface/atomic.py

@@ -237,7 +237,7 @@
     "Am": 1.8,
     "Cm": 1.69,
 }
-METALS = [3, 4, 11, 12, 13, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
+METALS = [3, 4, 5, 11, 12, 13, 14, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
           37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 55, 56, 57,
           58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74,
           75, 76, 77, 78, 79, 80, 81, 82, 83, 87, 88, 89, 90, 91, 92, 93, 94,

lammps_interface/lammps_main.py

@@ -436,8 +436,13 @@ def split_graph(self):
 
     def assign_force_fields(self):
 
-        attr = {'graph':self.graph, 'cutoff':self.options.cutoff, 'h_bonding':self.options.h_bonding,
-                'keep_metal_geometry':self.options.fix_metal, 'bondtype':self.options.dreid_bond_type}
+        attr = {
+            'graph': self.graph,
+            'cutoff': self.options.cutoff,
+            'h_bonding': self.options.h_bonding,
+            'keep_metal_geometry': self.options.fix_metal,
+            'bondtype': self.options.dreid_bond_type}
+
         param = getattr(ForceFields, self.options.force_field)(**attr)
 
         self.special_commands += param.special_commands()
@@ -517,11 +522,12 @@ def compute_simulation_size(self):
                     print("Use <ixjxk> format")
                     print("Exiting...")
                     sys.exit()
-        self.supercell=supercell
+
+        self.supercell = supercell
         if np.any(np.array(supercell) > 1):
             print("Re-sizing to a %i x %i x %i supercell. "%(supercell))
 
-            #TODO(pboyd): apply to subgraphs as well, if requested.
+            # TO DO (pboyd): apply to subgraphs as well, if requested.
             self.graph.build_supercell(supercell, self.cell)
             molcount = 0
             if self.subgraphs:
@@ -1140,8 +1146,11 @@ def construct_input_file(self):
                                      " ".join([self.unique_atom_types[key][1]['element']
                                                 for key in sorted(self.unique_atom_types.keys())])))
         elif self.options.dump_lammpstrj:
-            inp_str += "%-15s %s\n"%("dump","%s_lammpstrj all atom %i %s_mov.lammpstrj"%
-                                (self.name, self.options.dump_lammpstrj, self.name))
+            inp_str += "%-15s %s\n"%("dump","%s_lammpstrj all custom %i %s_mov.lammpstrj id element xs ys zs"% (self.name, self.options.dump_lammpstrj, self.name))
+
+            elements = " ".join([self.unique_atom_types[key][1]['element'] for key in sorted(self.unique_atom_types.keys())])
+            inp_str += "%-15s %s\n"%("dump_modify", "%s_lammpstrj sort id"% (self.name))
+            inp_str += "%-15s %s\n"%("dump_modify", "%s_lammpstrj element %s"% (self.name, elements))
 
             # in the meantime we need to map atom id to element that will allow us to
             # post-process the lammpstrj file and create a cif out of each
@@ -1152,52 +1161,71 @@ def construct_input_file(self):
             f.close()
 
         if (self.options.minimize):
-            box_min = "aniso"
+            # This part has been heavily modified by Felipe Lopes de Oliveira (FLO)
+            box_min = "tri"
             min_style = "cg"
-            min_eval = 1e-6   # HKUST-1 will not minimize past 1e-11
-            max_iterations = 100000 # if the minimizer can't reach a minimum in this many steps,
-                                    # change the min_eval to something higher.
-            #inp_str += "%-15s %s\n"%("min_style","fire")
-            #inp_str += "%-15s %i %s\n"%("compute", 1, "all msd com yes")
-            #inp_str += "%-15s %-10s %s\n"%("variable", "Dx", "equal c_1[1]")
-            #inp_str += "%-15s %-10s %s\n"%("variable", "Dy", "equal c_1[2]")
-            #inp_str += "%-15s %-10s %s\n"%("variable", "Dz", "equal c_1[3]")
-            #inp_str += "%-15s %-10s %s\n"%("variable", "MSD", "equal c_1[4]")
-            #inp_str += "%-15s %s %s\n"%("fix", "output all print 1", "\"$(vol),$(cella),$(cellb),$(cellc),${Dx},${Dy},${Dz},${MSD}\"" +
-            #                                " file %s.min.csv title \"Vol,CellA,CellB,CellC,Dx,Dy,Dz,MSD\" screen no"%(self.name))
-            inp_str += "%-15s %s\n"%("min_style", min_style)
-            inp_str += "%-15s %s\n"%("print", "\"MinStep,CellMinStep,AtomMinStep,FinalStep,Energy,EDiff\"" +
-                                              " file %s.min.csv screen no"%(self.name))
+            min_eval = self.options.min_eval   # HKUST-1 will not minimize past 1e-11
+            max_iterations = 100000  # if the minimizer can't reach a minimum in this many steps, change the min_eval to something higher.
+
+            inp_str += "#### Start of optimization ####\n\n"
+            inp_str += "%-15s %-10s %s\n"%("variable", "natoms", "equal count(all)")
+
+            inp_str += "%-15s %s\n"%("print", "\"Iter,PeTotal,DiffPe,PePerAtom,DiffPerAtom,Vol,F_max,F_Norm\"" +
+                                              " file %s.min.csv screen no\n"%(self.name))
+
+            inp_str += "# Loop Initiation\n\n"
+
             inp_str += "%-15s %-10s %s\n"%("variable", "min_eval", "equal %.2e"%(min_eval))
             inp_str += "%-15s %-10s %s\n"%("variable", "prev_E", "equal %.2f"%(50000.)) # set unreasonably high for first loop
+            inp_str += "%-15s %-10s %s\n"%("variable", "prev_E_atom", "equal %.2f"%(50000.)) # set unreasonably high for first loop
             inp_str += "%-15s %-10s %s\n"%("variable", "iter", "loop %i"%(max_iterations))
+
+            inp_str += "%-15s %-10s %s\n"%("variable", "F_max_comp", "equal fmax")
+            inp_str += "%-15s %-10s %s\n"%("variable", "F_norm_max", "equal fnorm")
+
             inp_str += "%-15s %s\n"%("label", "loop_min")
 
             fix = self.fixcount()
-            inp_str += "%-15s %s\n"%("min_style", min_style)
-            inp_str += "%-15s %s\n"%("fix","%i all box/relax %s 0.0 vmax 0.01"%(fix, box_min))
-            inp_str += "%-15s %s\n"%("minimize","1.0e-15 1.0e-15 10000 100000")
-            inp_str += "%-15s %s\n"%("unfix", "%i"%fix)
-            inp_str += "%-15s %s\n"%("min_style", "fire")
-            inp_str += "%-15s %-10s %s\n"%("variable", "tempstp", "equal $(step)")
-            inp_str += "%-15s %-10s %s\n"%("variable", "CellMinStep", "equal ${tempstp}")
-            inp_str += "%-15s %s\n"%("minimize","1.0e-15 1.0e-15 10000 100000")
-            inp_str += "%-15s %-10s %s\n"%("variable", "AtomMinStep", "equal $(step)")
-            inp_str += "%-15s %-10s %s\n"%("variable", "temppe", "equal $(pe)")
-            inp_str += "%-15s %-10s %s\n"%("variable", "min_E", "equal abs(${prev_E}-${temppe})")
-            inp_str += "%-15s %s\n"%("print", "\"${iter},${CellMinStep},${AtomMinStep},${AtomMinStep}," +
-                                              "$(pe),${min_E}\"" +
+
+            inp_str += "\n    # --- Stage 1: Box Relaxation ---\n"
+
+            inp_str += "    %-15s %s\n"%("min_style", min_style)
+            inp_str += "    %-15s\n"%("min_modify dmax 0.01")
+            inp_str += "    %-15s %s\n"%("fix","%i all box/relax %s 0.0 vmax 0.01"%(fix, box_min))
+            inp_str += "    %-15s %s\n"%("minimize","%.2e %.2e 5000 10000"%(min_eval, min_eval))
+            inp_str += "    %-15s %s\n"%("unfix", "%i"%fix)
+
+            inp_str += "\n    # --- Stage 2: Atom Relaxation ---\n"
+            inp_str += "    %-15s %s\n"%("min_style", "fire")
+            inp_str += "    %-15s\n"%("min_modify dmax 0.01")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "tempstp", "equal $(step)")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "CellMinStep", "equal ${tempstp}")
+            inp_str += "    %-15s %s\n"%("minimize","%.2e %.2e 10000 20000"%(min_eval, min_eval))
+
+            inp_str += "\n    # --- Convergence Check (Per Atom) ---\n"
+            inp_str += "    %-15s %-10s %s\n"%("variable", "curr_pe", "equal pe")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "curr_pe_atom", "equal v_curr_pe/v_natoms")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "vol", "equal vol")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "diff_total", "equal abs(${prev_E}-${curr_pe})")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "diff_atom", "equal abs(${prev_E_atom}-${curr_pe_atom})")
+
+            inp_str += "    %-15s %s\n\n"%("print", "\"${iter},${curr_pe},${diff_total},${curr_pe_atom},${diff_atom},${vol},${F_max_comp},${F_norm_max}\"" +
                                               " append %s.min.csv screen no"%(self.name))
 
-            inp_str += "%-15s %s\n"%("if","\"${min_E} < ${min_eval}\" then \"jump SELF break_min\"")
-            inp_str += "%-15s %-10s %s\n"%("variable", "prev_E", "equal ${temppe}")
-            inp_str += "%-15s %s\n"%("next", "iter")
+            inp_str += "    %-15s %-10s %s\n"%("variable", "prev_E", "equal ${curr_pe}")
+            inp_str += "    %-15s %-10s %s\n\n"%("variable", "prev_E_atom", "equal ${curr_pe_atom}")
+
+            inp_str += "    %-15s %s\n"%("if","\"${diff_atom} < ${min_eval}\" then \"jump SELF break_min\"")
+
+            inp_str += "    %-15s %s\n"%("next", "iter")
             inp_str += "%-15s %s\n"%("jump", "SELF loop_min")
             inp_str += "%-15s %s\n"%("label", "break_min")
+            inp_str += "\n#### End of optimization ####\n\n"
+
+            # Write the final relaxed structure to a file
+            inp_str += "%-15s %s\n"%("write_data", "data.optimized")
+
 
-           # inp_str += "%-15s %s\n"%("unfix", "output")
-        # delete bond types etc, for molecules that are rigid
-
         if (self.options.thermal_anneal):
             box_min = "aniso"
             min_style = "cg"
@@ -1244,29 +1272,6 @@ def construct_input_file(self):
                                                 for key in sorted(self.unique_atom_types.keys())])))
             inp_str += "run 0\n"
 
-        if (self.options.relax):
-            box_min = "aniso"
-            min_style = "cg"
-            inp_str += "run 0\n"
-            inp_str += "%-15s %-10s %s\n"%("variable", "inputpe", "equal $(pe)")
-            inp_str += "%-15s %s\n"%("min_style", min_style)
-            fix = self.fixcount()
-            inp_str += "%-15s %-10s %s\n"%("variable", "tempstp", "equal $(step)")
-            inp_str += "%-15s %s\n"%("minimize","1.0e-15 1.0e-15 10000 100000")
-            inp_str += "%-15s %-10s %s\n"%("variable", "AtomMinStep", "equal $(step)")
-            inp_str += "%-15s %-10s %s\n"%("variable", "relaxedpe", "equal $(pe)")
-            inp_str += "%-15s %-10s %s\n"%("variable", "min_E", "equal abs(${inputpe}-${relaxedpe})")
-            inp_str += "%-15s %s\n"%("print", "\"${AtomMinStep}," +
-                                              "$(pe),${min_E}\"" +
-                                              " append %s.min.csv screen no"%(self.name))
-            inp_str += "%-15s %s\n"%("dump","%s_relax all custom 1 relaxed_%s.lammpstrj element xs ys zs"%
-                            (self.name, self.name))
-            inp_str += "%-15s %s\n"%("dump_modify", "%s_relax element %s"%(
-                                     self.name,
-                                     " ".join([self.unique_atom_types[key][1]['element']
-                                                for key in sorted(self.unique_atom_types.keys())])))
-            inp_str += "run 0\n"
-
         for mol in sorted(self.molecule_types.keys()):
             rep = self.subgraphs[self.molecule_types[mol][0]]
             if rep.rigid: