Skip to content

Commit d6c52a8

Browse files
author
mlund
committed
xtcfile from cli
1 parent 847e30f commit d6c52a8

File tree

2 files changed

+44
-34
lines changed

2 files changed

+44
-34
lines changed

src/icoscan.rs

Lines changed: 32 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,7 @@ pub fn do_icoscan(
6868
temperature: &f64,
6969
pmf_file: &PathBuf,
7070
disktable: &Option<PathBuf>,
71+
xtcfile: &Option<PathBuf>,
7172
) -> std::result::Result<(), anyhow::Error> {
7273
let table = Table6D::from_resolution(rmin, rmax, dr, angle_resolution)?;
7374
let n_vertices = table.get(rmin)?.get(0.0)?.len();
@@ -118,34 +119,36 @@ pub fn do_icoscan(
118119
});
119120

120121
// Write oriented structures to trajectory file
121-
log::info!("Writing trajectory file traj.xtc");
122-
let mut traj = xdrfile::XTCTrajectory::open_write("traj.xtc")?;
123-
let mut frame_cnt: usize = 0;
124-
let n = r_and_omega.len();
125-
r_and_omega
126-
.into_iter()
127-
.progress_count(n as u64)
128-
.for_each(|(r, omega)| {
129-
table
130-
.get_icospheres(r, omega) // remaining 4D
131-
.expect("invalid (r, omega) value")
132-
.flat_iter()
133-
.for_each(|(pos_a, pos_b, _data_b)| {
134-
let (oriented_a, oriented_b) =
135-
orient_structures(r, omega, *pos_a, *pos_b, &ref_a, &ref_b);
136-
let mut frame = xdrfile::Frame::new();
137-
frame.step = frame_cnt;
138-
frame.time = frame_cnt as f32;
139-
frame_cnt += 1;
140-
frame.coords = oriented_a
141-
.pos
142-
.iter()
143-
.chain(oriented_b.pos.iter())
144-
.map(|&p| [p.x as f32, p.y as f32, p.z as f32])
145-
.collect();
146-
traj.write(&frame).expect("Failed to write XTC frame");
147-
});
148-
});
122+
if let Some(xtcfile) = xtcfile {
123+
log::info!("Writing trajectory file {}", xtcfile.display());
124+
let mut traj = xdrfile::XTCTrajectory::open_write(xtcfile)?;
125+
let mut frame_cnt: usize = 0;
126+
let n = r_and_omega.len();
127+
r_and_omega
128+
.into_iter()
129+
.progress_count(n as u64)
130+
.for_each(|(r, omega)| {
131+
table
132+
.get_icospheres(r, omega) // remaining 4D
133+
.expect("invalid (r, omega) value")
134+
.flat_iter()
135+
.for_each(|(pos_a, pos_b, _data_b)| {
136+
let (oriented_a, oriented_b) =
137+
orient_structures(r, omega, *pos_a, *pos_b, &ref_a, &ref_b);
138+
let mut frame = xdrfile::Frame::new();
139+
frame.step = frame_cnt;
140+
frame.time = frame_cnt as f32;
141+
frame_cnt += 1;
142+
frame.coords = oriented_a
143+
.pos
144+
.iter()
145+
.chain(oriented_b.pos.iter())
146+
.map(|&p| [p.x as f32, p.y as f32, p.z as f32])
147+
.collect();
148+
traj.write(&frame).expect("Failed to write XTC frame");
149+
});
150+
});
151+
}
149152

150153
// Partition function contribution for single (r, omega) point
151154
// i.e. averaged over 4D angular space
@@ -163,7 +166,7 @@ pub fn do_icoscan(
163166
if let Some(savetable) = disktable {
164167
log::info!("Saving 6D table to {}", savetable.display());
165168
let mut stream = faunus::aux::open_compressed(savetable)?;
166-
for r in distances.iter() {
169+
for r in distances.iter().progress_count(distances.len() as u64) {
167170
table.stream_angular_space(*r, &mut stream)?;
168171
}
169172
}

src/main.rs

Lines changed: 12 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -125,6 +125,9 @@ enum Commands {
125125
/// Save table to disk (use .gz suffix for compression)
126126
#[arg(long)]
127127
savetable: Option<PathBuf>,
128+
/// Export XTC file with all poses
129+
#[arg(long)]
130+
xtcfile: Option<PathBuf>,
128131
},
129132
}
130133

@@ -144,6 +147,7 @@ fn do_scan(cmd: &Commands) -> Result<()> {
144147
fixed_dielectric,
145148
pmf_file,
146149
savetable,
150+
xtcfile,
147151
} = cmd
148152
else {
149153
bail!("Unknown command");
@@ -175,11 +179,13 @@ fn do_scan(cmd: &Commands) -> Result<()> {
175179
);
176180
let ref_a = Structure::from_xyz(mol1, topology.atomkinds())?;
177181
let ref_b = Structure::from_xyz(mol2, topology.atomkinds())?;
178-
let mut xyz_file = File::create("confout.xyz")?;
179-
ref_a
180-
.clone()
181-
.add(ref_b.clone())
182-
.to_xyz(&mut xyz_file, topology.atomkinds())?;
182+
if xtcfile.is_some() {
183+
log::info!("Exporting merged XYZ file with both initial structures: confout.xyz");
184+
ref_a
185+
.clone()
186+
.add(ref_b.clone())
187+
.to_xyz(&mut File::create("confout.xyz")?, topology.atomkinds())?;
188+
}
183189

184190
info!("{}", medium);
185191
info!(
@@ -216,6 +222,7 @@ fn do_scan(cmd: &Commands) -> Result<()> {
216222
temperature,
217223
pmf_file,
218224
savetable,
225+
xtcfile,
219226
)
220227
}
221228

0 commit comments

Comments
 (0)