Skip to content

Commit eb10b40

Browse files
srmnitcsrmnitc
committed
wip
1 parent 2bc5627 commit eb10b40

File tree

4 files changed

+800
-37
lines changed

4 files changed

+800
-37
lines changed

‎calphy/phase.py‎

Lines changed: 82 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1059,32 +1059,32 @@ def _reversible_scaling_forward(self, iteration=1):
10591059
l_values=self.calc.structural_monitoring.l_values,
10601060
cutoff=self.calc.structural_monitoring.cutoff,
10611061
detector_type="cumulative",
1062-
threshold=3,
1062+
threshold=3 / 2,
10631063
)
10641064

10651065
# Set up dump for equilibration trajectory to pre-train
1066-
pretrain_dump, pretrain_dump_id = sm.setup_pretraining_dump(
1067-
lmp=lmp, simfolder=self.simfolder, dump_interval=100
1068-
)
1069-
self.logger.info(
1070-
"Dumping equilibration trajectory every 100 steps for pre-training"
1071-
)
1066+
# pretrain_dump, pretrain_dump_id = sm.setup_pretraining_dump(
1067+
# lmp=lmp, simfolder=self.simfolder, dump_interval=100
1068+
# )
1069+
# self.logger.info(
1070+
# "Dumping equilibration trajectory every 100 steps for pre-training"
1071+
# )
10721072

10731073
self.logger.info(f"Starting equilibration with constrained com: {iteration}")
10741074
lmp.command("run %d" % self.calc.n_equilibration_steps)
10751075
self.logger.info(f"Finished equilibration with constrained com: {iteration}")
10761076

10771077
# Pre-train FrameAccumulator on equilibration trajectory
1078-
if accumulator is not None and pretrain_dump is not None:
1079-
lmp.command(f"undump {pretrain_dump_id}")
1080-
1081-
success = sm.pretrain_accumulator(
1082-
accumulator=accumulator, pretrain_dump=pretrain_dump, logger=self.logger
1083-
)
1078+
# if accumulator is not None and pretrain_dump is not None:
1079+
# lmp.command(f"undump {pretrain_dump_id}")
1080+
#
1081+
# success = sm.pretrain_accumulator(
1082+
# accumulator=accumulator, pretrain_dump=pretrain_dump, logger=self.logger
1083+
# )
10841084

1085-
# Clean up trajectory file
1086-
if os.path.exists(pretrain_dump):
1087-
os.remove(pretrain_dump)
1085+
# # Clean up trajectory file
1086+
# if os.path.exists(pretrain_dump):
1087+
# os.remove(pretrain_dump)
10881088

10891089
lmp.command("variable flambda equal ramp(${li},${lf})")
10901090
lmp.command("variable blambda equal ramp(${lf},${li})")
@@ -1222,10 +1222,6 @@ def _reversible_scaling_forward(self, iteration=1):
12221222
self.logger.info(
12231223
f"Temperature-based monitoring: {len(checkpoint_data)-1} blocks"
12241224
)
1225-
for i, cp in enumerate(checkpoint_data):
1226-
self.logger.info(
1227-
f" Checkpoint {i}: T={cp['temp']:.2f}K, λ={cp['lambda']:.4f}, step={cp['step']}"
1228-
)
12291225

12301226
# Open file to save monitoring data
12311227
monitor_file = os.path.join(
@@ -1266,15 +1262,76 @@ def _reversible_scaling_forward(self, iteration=1):
12661262
try:
12671263
atoms = sm.extract_atoms_from_lammps(lmp)
12681264

1265+
# here lets start a new LAMMPS instance and then relax it
1266+
# create lammps object
1267+
# Get detector statistics for logging
1268+
lambda_value = li + (lf - li) * current_step / total_steps
1269+
apparent_temp = t0 / lambda_value
1270+
1271+
self.logger.info(
1272+
f" Checkpoint {block_idx}/{len(checkpoint_data)-1}: T={apparent_temp:.2f}K, λ={lambda_value:.4f}"
1273+
)
1274+
self.logger.info(
1275+
"Starting quick relaxation for structural monitoring at block %d"
1276+
% block_idx
1277+
)
1278+
lmpmin = ph.create_object(
1279+
self.cores,
1280+
self.simfolder,
1281+
self.calc.md.timestep,
1282+
self.calc.md.cmdargs,
1283+
self.calc.md.init_commands,
1284+
)
1285+
1286+
lmpmin.command("echo log")
1287+
lmpmin.command("variable li equal %f" % li)
1288+
lmpmin.command("variable lf equal %f" % lf)
1289+
1290+
lmpmin.command(
1291+
f"pair_style {self.calc._pair_style_with_options[0]}"
1292+
)
1293+
1294+
# read in conf file
1295+
# conf = os.path.join(self.simfolder, "conf.equilibration.dump")
1296+
from ase.io import write
1297+
1298+
write(
1299+
os.path.join(self.simfolder, "temp_monitoring_frame.data"),
1300+
atoms,
1301+
format="lammps-data",
1302+
atom_style="atomic",
1303+
)
1304+
conf = os.path.join(self.simfolder, "temp_monitoring_frame.data")
1305+
lmpmin = ph.read_data(lmpmin, conf)
1306+
1307+
# set up potential
1308+
lmpmin.command(f"pair_coeff {self.calc.pair_coeff[0]}")
1309+
lmpmin = ph.set_mass(lmpmin, self.calc)
1310+
1311+
# now we do a quick relaxation to get the structure back to a reasonable state before calculating the stats
1312+
lmpmin.command("fix ensemble all box/relax aniso 0.0")
1313+
lmpmin.command(f"minimize 1e-5 1e-5 100000 100000")
1314+
lmpmin.command("run 0")
1315+
1316+
atoms = sm.extract_atoms_from_lammps(lmpmin)
1317+
lmpmin.close()
1318+
self.logger.info(
1319+
"Completed quick relaxation for structural monitoring at block %d"
1320+
% block_idx
1321+
)
1322+
write(
1323+
os.path.join(
1324+
self.simfolder, f"temp_monitoring_frame_{block_idx}.data"
1325+
),
1326+
atoms,
1327+
format="lammps-data",
1328+
atom_style="atomic",
1329+
)
12691330
# Process frame and get statistics
12701331
stats = sm.process_monitoring_frame(
12711332
accumulator=accumulator, atoms=atoms, logger=self.logger
12721333
)
12731334

1274-
# Get detector statistics for logging
1275-
lambda_value = li + (lf - li) * current_step / total_steps
1276-
apparent_temp = t0 / lambda_value
1277-
12781335
# Write to monitoring file
12791336
monitor_out.write(
12801337
f"{block_idx} {current_step} {lambda_value:.6f} {apparent_temp:.2f} "
@@ -1497,7 +1554,7 @@ def _reversible_scaling_backward(self, iteration=1):
14971554
l_values=self.calc.structural_monitoring.l_values,
14981555
cutoff=self.calc.structural_monitoring.cutoff,
14991556
detector_type="cumulative",
1500-
threshold=3,
1557+
threshold=3 / 2,
15011558
)
15021559

15031560
# Set up dump for equilibration trajectory to pre-train

‎calphy/structural_monitoring.py‎

Lines changed: 22 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,8 @@ def process_monitoring_frame(accumulator, atoms, logger=None):
223223
# Get the latest distance
224224
if len(distances_arr) > 0:
225225
stats["distance"] = distances_arr[-1]
226+
else:
227+
stats["distance"] = np.nan
226228

227229
# Check if current frame is flagged
228230
current_frame_idx = len(frames_arr) - 1 if len(frames_arr) > 0 else -1
@@ -249,16 +251,20 @@ def process_monitoring_frame(accumulator, atoms, logger=None):
249251
"std": float(np.std(finite)),
250252
}
251253

252-
stats["mean"] = (
253-
detector_stats.get("mean", -1.0)
254-
if detector_stats.get("mean") is not None
255-
else -1.0
256-
)
257-
stats["std"] = (
258-
detector_stats.get("std", -1.0)
259-
if detector_stats.get("std") is not None
260-
else -1.0
261-
)
254+
# If no previous data, use current value for mean, 0 for std
255+
if detector_stats.get("mean") is not None:
256+
stats["mean"] = detector_stats["mean"]
257+
elif stats["distance"] != -1.0:
258+
stats["mean"] = stats["distance"]
259+
else:
260+
stats["mean"] = np.nan
261+
262+
if detector_stats.get("std") is not None:
263+
stats["std"] = detector_stats["std"]
264+
elif stats["distance"] != -1.0:
265+
stats["std"] = 0.0
266+
else:
267+
stats["std"] = np.nan
262268

263269
except Exception as e:
264270
if logger:
@@ -311,8 +317,12 @@ def generate_monitoring_plot(monitor_file, output_file, iteration=1):
311317
mean_arr = data[:, 6]
312318
std_arr = data[:, 7]
313319
else:
314-
mean_arr = np.array([np.mean(distance_arr[: i + 1]) for i in range(len(distance_arr))])
315-
std_arr = np.array([np.std(distance_arr[: i + 1]) for i in range(len(distance_arr))])
320+
mean_arr = np.array(
321+
[np.mean(distance_arr[: i + 1]) for i in range(len(distance_arr))]
322+
)
323+
std_arr = np.array(
324+
[np.std(distance_arr[: i + 1]) for i in range(len(distance_arr))]
325+
)
316326

317327
# Create plot
318328
fig, ax = plt.subplots(figsize=(10, 6))

0 commit comments

Comments
 (0)