Skip to content

Read bond types from CIF file. #76

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lipelopesoliveira wants to merge 12 commits into
base: master
Choose a base branch
from
Open

Read bond types from CIF file. #76

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
f9cb05a
Dump the optimized structure after minimization of relax as a "data.o…
lipelopesoliveira Apr 22, 2025
519eec0
Add the `bond_types_from_cif` option to read the bond types from CIF …
lipelopesoliveira Apr 22, 2025
2916177
Set bond_types_from_cif to False when no bonds are reported in CIF file
lipelopesoliveira Apr 22, 2025
76832f7
Refactor from_CIF function parameters for clarity
lipelopesoliveira Apr 22, 2025
7e10c60
Finishes the implementation of `ase_from_CIF` method to read the cif …
lipelopesoliveira Apr 22, 2025
5b0427d
Fix typo in variable name for distance checking in MolecularGraph ini…
lipelopesoliveira Apr 22, 2025
d372f8e
Improve trajectory writing to include atom types
lipelopesoliveira Jan 13, 2026
c4cc23a
Refactor dump_modify command to improve readability of element mappin…
lipelopesoliveira Jan 13, 2026
b556817
Refactor optimization on minimize
lipelopesoliveira Jan 13, 2026
04bdc95
Enhance sulfur force field type assignment based on graph degree
lipelopesoliveira Jan 13, 2026
06e86db
Refactor minimize command parameters to use dynamic evaluation for be…
lipelopesoliveira Jan 13, 2026
b463ff2
Add minimal energy change parameter for optimization convergence
lipelopesoliveira Jan 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 14 additions & 2 deletions lammps_interface/ForceFields.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2499,8 +2499,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']
Expand Down Expand Up @@ -3616,6 +3620,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"
Expand Down
8 changes: 8 additions & 0 deletions lammps_interface/InputHandler.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
129 changes: 66 additions & 63 deletions lammps_interface/lammps_main.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand All @@ -1152,52 +1161,69 @@ 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 %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 %-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"
Expand Down Expand Up @@ -1244,29 +1270,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:
Expand Down
Loading