Implement pair_potential function for pets and pcsaft (#24)

g-bauer · web-flow · commit 806a8e7558de · 2022-07-07T15:40:01.000+02:00
* Implement pair_potential function for pets and pcsaft

* fixes
diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs
@@ -166,6 +166,11 @@ impl FluidParameters for PcSaftFunctional {
 
 impl PairPotential for PcSaftFunctional {
     fn pair_potential(&self, i: usize, r: &Array1<f64>, _: f64) -> Array2<f64> {
-        unimplemented!()
+        let sigma_ij = &self.parameters.sigma_ij;
+        let eps_ij_4 = 4.0 * &self.parameters.epsilon_k_ij;
+        Array::from_shape_fn((self.parameters.m.len(), r.len()), |(j, k)| {
+            let att = (sigma_ij[[i, j]] / r[k]).powi(6);
+            eps_ij_4[[i, j]] * att * (att - 1.0)
+        })
     }
 }
diff --git a/src/pets/dft/mod.rs b/src/pets/dft/mod.rs
@@ -143,6 +143,16 @@ impl FluidParameters for PetsFunctional {
 
 impl PairPotential for PetsFunctional {
     fn pair_potential(&self, i: usize, r: &Array1<f64>, _: f64) -> Array2<f64> {
-        todo!()
+        let eps_ij_4 = 4.0 * self.parameters.epsilon_k_ij.clone();
+        let shift_ij = &eps_ij_4 * (2.5.powi(-12) - 2.5.powi(-6));
+        let rc_ij = 2.5 * &self.parameters.sigma_ij;
+        Array::from_shape_fn((self.parameters.sigma.len(), r.len()), |(j, k)| {
+            if r[k] > rc_ij[[i, j]] {
+                0.0
+            } else {
+                let att = (self.parameters.sigma_ij[[i, j]] / r[k]).powi(6);
+                eps_ij_4[[i, j]] * att * (att - 1.0) - shift_ij[[i, j]]
+            }
+        })
     }
 }

Original file line number	Diff line number	Diff line change
`@@ -166,6 +166,11 @@ impl FluidParameters for PcSaftFunctional {`
`166`	`166`
`167`	`167`	`impl PairPotential for PcSaftFunctional {`
`168`	`168`	`fn pair_potential(&self, i: usize, r: &Array1<f64>, _: f64) -> Array2<f64> {`
`169`		`- unimplemented!()`
	`169`	`+ let sigma_ij = &self.parameters.sigma_ij;`
	`170`	`+ let eps_ij_4 = 4.0 * &self.parameters.epsilon_k_ij;`
	`171`	`+ Array::from_shape_fn((self.parameters.m.len(), r.len()), \|(j, k)\| {`
	`172`	`+ let att = (sigma_ij[[i, j]] / r[k]).powi(6);`
	`173`	`+ eps_ij_4[[i, j]] * att * (att - 1.0)`
	`174`	`+ })`
`170`	`175`	`}`
`171`	`176`	`}`