Merge pull request #19 from benmagruder/bugfix

Fixed bug that gives pressure=NaN when a species has phi=0
dmorse · May 4, 2022 · c01e82f · c01e82f
2 parents b9409d7 + 48d9f7f
commit c01e82f
Showing 1 changed file with 6 additions and 2 deletions.
diff --git 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