From 48d9f7fb6c4885c0506d0287503cf575ab723cee Mon Sep 17 00:00:00 2001 From: Ben Magruder Date: Tue, 3 May 2022 14:55:18 -0500 Subject: [PATCH] Fixed bug that gives pressure=nan when a species has phi=0 --- src/scf_mod.fp.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/scf_mod.fp.f90 b/src/scf_mod.fp.f90 index 362f232..b61f187 100644 --- a/src/scf_mod.fp.f90 +++ b/src/scf_mod.fp.f90 @@ -872,10 +872,14 @@ subroutine free_energy(N, rho, omega, phi_chain, mu_chain, & if (present(pressure)) then pressure = -f_Helmholtz do i = 1, N_chain - pressure = pressure + mu_chain(i)*phi_chain(i)/chain_length(i) + if ( phi_chain(i) > 1.0E-8) then + pressure = pressure + mu_chain(i)*phi_chain(i)/chain_length(i) + end if end do do i = 1, N_solvent - pressure = pressure + mu_solvent(i)*phi_solvent(i)/solvent_size(i) + if ( phi_solvent(i) > 1.0E-8) then + pressure = pressure + mu_solvent(i)*phi_solvent(i)/solvent_size(i) + end if end do end if