Skip to content

Commit

Permalink
Added PMV density; full writer
Browse files Browse the repository at this point in the history 
PMV density needs a choice between scalar or density-mapped isothem-comp
Everything should be composed into a calculator struct now
  • Loading branch information
@2AUK
2AUK committed Nov 7, 2023
1 parent c06596e commit e63053f
Show file tree
Hide file tree
Showing 5 changed files with 293 additions and 101 deletions.
15 changes: 15 additions & 0 deletions SPCE_DRISM_methane_UA.td
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
Thermodynamics:
Temperature: 283.15 K
Total Density: 0.03333 (1/A^3)
Isothermal Compressibility: 1.993345317372601
Pressure: 0.23562577083196018
Molecular KB PMV: -71.8164512050989 A^3
-43.24786691571056 cm^3/mol
RISM KB PMV: 52.08188812018311 A^3
31.363713025974267 cm^3/mol
Solvation Free Energies (kcal/mol):
HNC: 9.171794554957005
KH: 8.015539657449239
GF: 5.4701122255142165
PW: 5.544899105410485
PC+: 12.271835034702054
2 changes: 1 addition & 1 deletion pyrism/plotter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
def plot():
plt.rcParams["text.usetex"] = False
plt.rcParams["font.size"] = 14
df = pd.read_csv(sys.argv[1], sep=",", skiprows=[0])
df = pd.read_csv(sys.argv[1], sep=",")

r = df.iloc[:, 0].to_numpy()
grs = []
Expand Down
237 changes: 173 additions & 64 deletions src/io/writer.rs
View file
Original file line number Diff line number Diff line change
@@ -1,111 +1,220 @@
use crate::data::solution::Solutions;
use crate::grids::radial_grid::Grid;
use crate::{data::solution::Solutions, thermodynamics::thermo::Thermodynamics};
use csv::{QuoteStyle, WriterBuilder};
use ndarray::Array3;
use itertools::Itertools;
use ndarray::{Array1, Array3};
use std::fs::File;
use std::io::prelude::*;
use std::path::Path;

pub struct RISMWriter<'a> {
pub name: String,
pub data: &'a Solutions,
pub thermo: &'a Thermodynamics,
}

impl<'a> RISMWriter<'a> {
pub fn new(name: &String, data: &'a Solutions) -> Self {
pub fn new(name: &String, data: &'a Solutions, thermo: &'a Thermodynamics) -> Self {
RISMWriter {
name: name.clone(),
data,
thermo,
}
}

pub fn write(&self) -> std::io::Result<()> {
let vv = &self.data.vv;
let tr = &vv.correlations.tr;
let cr = &vv.correlations.cr;
let gr = &vv.correlations.tr + &vv.correlations.cr + 1.0;
self.write_file_vv(&gr, "gvv".to_string())?;
self.write_file_vv(&vv.correlations.tr, "tvv".to_string())?;
self.write_file_vv(&vv.correlations.cr, "cvv".to_string())?;
// Set up grid for dumping data with
let grid = Grid::new(vv.data_config.npts, vv.data_config.radius);
// Construct header row for .csv file from atom-site names
let mut vv_headers: Vec<String> = vv
.data_config
.solvent_atoms
.iter()
.cartesian_product(vv.data_config.solvent_atoms.iter())
.map(|(iat, jat)| format!("{}-{}", &iat.atom_type, &jat.atom_type))
.collect();
vv_headers.insert(0, "r(A)".to_string());
// Write correlation functions
self.write_correlation(&gr, &grid.rgrid, &vv_headers, "gvv".to_string(), gr.dim())?;
self.write_correlation(cr, &grid.rgrid, &vv_headers, "cvv".to_string(), gr.dim())?;
self.write_correlation(tr, &grid.rgrid, &vv_headers, "tvv".to_string(), gr.dim())?;

// Check for a uv problem and write files if present
match self.data.uv {
Some(ref uv) => {
let gr = &uv.correlations.tr + &uv.correlations.cr + 1.0;
self.write_file_uv(&gr, "guv".to_string())?;
self.write_file_uv(&uv.correlations.tr, "tuv".to_string())?;
self.write_file_uv(&uv.correlations.cr, "cuv".to_string())?;
let tr = &uv.correlations.tr;
let cr = &uv.correlations.cr;
let gr = tr + cr + 1.0;
let mut uv_headers: Vec<String> = uv
.data_config
.solvent_atoms
.iter()
.cartesian_product(vv.data_config.solvent_atoms.iter())
.map(|(iat, jat)| format!("{}-{}", &iat.atom_type, &jat.atom_type))
.collect();

uv_headers.insert(0, "r(A)".to_string());
self.write_correlation(&gr, &grid.rgrid, &uv_headers, "gvv".to_string(), gr.dim())?;
self.write_correlation(cr, &grid.rgrid, &uv_headers, "cvv".to_string(), gr.dim())?;
self.write_correlation(tr, &grid.rgrid, &uv_headers, "tvv".to_string(), gr.dim())?;

//Write thermodynamics to file
let td_path = format!("{}.td", self.name);
let mut file = File::create(td_path)?;
file.write_all(self.thermo.to_string().as_bytes())?;

//Write densities to file
self.write_density(&grid.rgrid, gr.dim())?;
Ok(())
}
None => {
//Write solvent-solvent thermodynamics only
let td_path = format!("{}.td", self.name);
let mut file = File::create(td_path)?;
file.write_all(self.thermo.to_string().as_bytes())?;
Ok(())
}
None => Ok(()),
}
}

fn write_file_uv(&self, func: &Array3<f64>, ext: String) -> std::io::Result<()> {
let vv = &self.data.vv;
let uv = self.data.uv.as_ref().unwrap();
let grid = Grid::new(vv.data_config.npts, vv.data_config.radius);
let path_str = format!("{}_rust.{}", self.name, ext);
let path = Path::new(&path_str);
fn write_correlation(
&self,
func: &Array3<f64>,
grid: &Array1<f64>,
header: &[String],
ext: String,
_shape @ (npts, ns1, ns2): (usize, usize, usize),
) -> std::io::Result<()> {
let path_str = &(format!("{}.{}", self.name, ext)).to_string();
let path = Path::new(path_str);
let mut wtr = WriterBuilder::new()
.comment(Some(b'#'))
.flexible(true)
.quote_style(QuoteStyle::Never)
.from_path(path)?;
wtr.write_record(&["#input"])?;
let mut header_row = Vec::new();
header_row.push("r(A)".to_string());
// generate header row
for iat in uv.data_config.solvent_atoms.iter() {
for jat in vv.data_config.solvent_atoms.iter() {
let header_string = format!("{}-{}", &iat.atom_type, &jat.atom_type);
header_row.push(header_string);
}
}
wtr.write_record(header_row.as_slice())?;
// tabulate data
for i in 0..uv.data_config.npts {
let mut data_row = Vec::new();
data_row.push(grid.rgrid[[i]].to_string());
for j in 0..uv.data_config.nsu.unwrap() {
for k in 0..vv.data_config.nsv {
data_row.push(func[[i, j, k]].to_string());
wtr.write_record(header)?;
for i in 0..npts {
let mut data = Vec::new();
data.push(grid[[i]].to_string());
for j in 0..ns1 {
for k in 0..ns2 {
data.push(func[[i, j, k]].to_string());
}
}
wtr.write_record(&data_row[..])?;
wtr.write_record(data)?;
}
wtr.flush()?;
Ok(())
}

fn write_file_vv(&self, func: &Array3<f64>, ext: String) -> std::io::Result<()> {
let vv = &self.data.vv;
let grid = Grid::new(vv.data_config.npts, vv.data_config.radius);
let path_str = format!("{}_rust.{}", self.name, ext);
let path = Path::new(&path_str);
fn write_density(
&self,
grid: &Array1<f64>,
_shape @ (npts, _, _): (usize, usize, usize),
) -> std::io::Result<()> {
let path_str = &(format!("{}.duv", self.name)).to_string();
let path = Path::new(path_str);
let mut wtr = WriterBuilder::new()
.comment(Some(b'#'))
.flexible(true)
.quote_style(QuoteStyle::Never)
.from_path(path)?;
wtr.write_record(&["#input"])?;
let mut header_row = Vec::new();
header_row.push("r(A)".to_string());
// generate header row
for iat in vv.data_config.solvent_atoms.iter() {
for jat in vv.data_config.solvent_atoms.iter() {
let header_string = format!("{}-{}", &iat.atom_type, &jat.atom_type);
header_row.push(header_string);
}
let header = vec![
"r(A)".to_string(),
"HNC".to_string(),
"KH".to_string(),
"GF".to_string(),
"PW".to_string(),
"PMV".to_string(),
];
wtr.write_record(header.as_slice())?;
let densities = self.thermo.sfed.as_ref().unwrap();
for i in 0..npts {
let data = vec![
grid[[i]].to_string(),
densities.hypernettedchain[[i]].to_string(),
densities.kovalenko_hirata[[i]].to_string(),
densities.gaussian_fluctuations[[i]].to_string(),
densities.partial_wave[[i]].to_string(),
densities.partial_molar_volume[[i]].to_string(),
];
wtr.write_record(data.as_slice())?;
}
wtr.write_record(header_row.as_slice())?;
// tabulate data
for i in 0..vv.data_config.npts {
let mut data_row = Vec::new();
data_row.push(grid.rgrid[[i]].to_string());
for j in 0..vv.data_config.nsv {
for k in 0..vv.data_config.nsv {
data_row.push(func[[i, j, k]].to_string());
}
}
wtr.write_record(&data_row[..])?;
}
wtr.flush()?;
Ok(())
}

// fn write_file_uv(&self, func: &Array3<f64>, ext: String) -> std::io::Result<()> {
// let vv = &self.data.vv;
// let uv = self.data.uv.as_ref().unwrap();
// let path_str = format!("{}_rust.{}", self.name, ext);
// let path = Path::new(&path_str);
// let mut wtr = WriterBuilder::new()
// .comment(Some(b'#'))
// .flexible(true)
// .quote_style(QuoteStyle::Never)
// .from_path(path)?;
// wtr.write_record(&["#input"])?;
// let mut header_row = Vec::new();
// header_row.push("r(A)".to_string());
// // generate header row
// for iat in uv.data_config.solvent_atoms.iter() {
// for jat in vv.data_config.solvent_atoms.iter() {
// let header_string = format!("{}-{}", &iat.atom_type, &jat.atom_type);
// header_row.push(header_string);
// }
// }
// wtr.write_record(header_row.as_slice())?;
// // tabulate data
// for i in 0..uv.data_config.npts {
// let mut data_row = Vec::new();
// data_row.push(grid.rgrid[[i]].to_string());
// for j in 0..uv.data_config.nsu.unwrap() {
// for k in 0..vv.data_config.nsv {
// data_row.push(func[[i, j, k]].to_string());
// }
// }
// wtr.write_record(&data_row[..])?;
// }
// wtr.flush()?;
// Ok(())
// }
//
// fn write_file_vv(&self, func: &Array3<f64>, ext: String) -> std::io::Result<()> {
// let vv = &self.data.vv;
// let grid = Grid::new(vv.data_config.npts, vv.data_config.radius);
// let path_str = format!("{}_rust.{}", self.name, ext);
// let path = Path::new(&path_str);
// let mut wtr = WriterBuilder::new()
// .comment(Some(b'#'))
// .flexible(true)
// .quote_style(QuoteStyle::Never)
// .from_path(path)?;
// wtr.write_record(&["#input"])?;
// let mut header_row = Vec::new();
// header_row.push("r(A)".to_string());
// // generate header row
// for iat in vv.data_config.solvent_atoms.iter() {
// for jat in vv.data_config.solvent_atoms.iter() {
// let header_string = format!("{}-{}", &iat.atom_type, &jat.atom_type);
// header_row.push(header_string);
// }
// }
// wtr.write_record(header_row.as_slice())?;
// // tabulate data
// for i in 0..vv.data_config.npts {
// let mut data_row = Vec::new();
// data_row.push(grid.rgrid[[i]].to_string());
// for j in 0..vv.data_config.nsv {
// for k in 0..vv.data_config.nsv {
// data_row.push(func[[i, j, k]].to_string());
// }
// }
// wtr.write_record(&data_row[..])?;
// }
// wtr.flush()?;
// Ok(())
// }
}
21 changes: 11 additions & 10 deletions src/main.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,16 @@ fn main() -> Result<(), lexopt::Error> {
let args = parse_args()?;
let mut driver = RISMDriver::from_toml(&args.input_file);
let solutions = driver.execute(args.verbosity, args.compress);

let wv = driver.solvent.borrow().wk.clone();
let wu = {
match driver.solute {
Some(v) => Some(v.borrow().wk.clone()),
None => None,
}
};
let td = TDDriver::new(&solutions, wv, wu);
let thermo = td.execute();
let writer = RISMWriter::new(
&args
.input_file
Expand All @@ -67,17 +77,8 @@ fn main() -> Result<(), lexopt::Error> {
.unwrap()
.to_string(),
&solutions,
&thermo,
);
writer.write().unwrap();
let wv = driver.solvent.borrow().wk.clone();
let wu = {
match driver.solute {
Some(v) => Some(v.borrow().wk.clone()),
None => None,
}
};
let td = TDDriver::new(solutions, wv, wu);
let thermo = td.execute();
println!("{}", thermo);
Ok(())
}

0 comments on commit e63053f

Please sign in to comment.