1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
use std::iter;
use anyhow::Result;
use log::info;
use molar::prelude::*;
pub(crate) fn command_tip3_to_tip4(file: &str, outfile: &str) -> Result<()> {
info!("Loading file '{file}'...");
let inp = System::from_file(file)?;
// Make output system
let mut out = System::default();
// Select water
// We don't modify input, so all selections could be bound
let water = inp.select_bound("resname TIP3")?;
// For correct re-assembly of the system
// select what is before and what is after water
let w_first = water.first_index();
let w_last = water.last_index();
let sel_before = inp.select_bound(0..w_first)?;
let sel_after = inp.select_bound(w_last + 1..inp.len())?;
// Add before selection
out.append(&sel_before)?;
// Now go over water molecules one by one
for mol in water.split_resindex_bound() {
// TIP3 is arranged as O->H->H
// so atom 0 is O, atoms 1 and 2 are H
// Get cooridnates
let o_pos = *mol.get_pos(0).unwrap();
let h1_pos = *mol.get_pos(1).unwrap();
let h2_pos = *mol.get_pos(2).unwrap();
// Get center of masses of H
let hc = 0.5 * (h1_pos.coords + h2_pos.coords);
// Unit vector from o to hc
let v = (hc - o_pos.coords).normalize();
// Position of the M dummy particle in TIP4
let m_pos = o_pos + v * 0.01546;
// Dummy atom M
let m_at = Atom {
name: AtomStr::try_from_str("M").unwrap(),
..mol.first_particle().atom.clone()
};
println!("{:?}", m_at);
// Add new converted water molecule
// We assume that the dummy is the last atom.
let added = out.append_atoms_coords(
mol.iter_atoms().chain(iter::once(&m_at)),
mol.iter_pos().chain(iter::once(&m_pos)),
)?;
// Change resname for added atoms
out.select_bound_mut(added)?.set_same_resname("TIP4");
}
// Add after selection
out.append(&sel_after)?;
// Transfer box
out.set_box_from(&inp);
// Write out new system
out.save(outfile)?;
Ok(())
}