!** array for fugacities
Allocate(gcmcparams%fuglist(nsims, nspc), STAT=error)
If (error/=0) Call allocErrDisplay(__FILE__,__LINE__,'fuglist')
!** Get the fugacity
Do i=1, nsims
!** Set the pressure of each component
Do j=1, gcmcparams%nspc
spc = gcmcparams%gcmcspc(j)%spc
!** pp is in kPa
pp = gcmcparams%gcmcspc(j)%fuglist(i)%pressure
! Write(*,*) spc, pp
Call eos_setpressure(gcmcparams%eosparams, spc, pp)
End Do
!** Get the fugacity of each component
Do j=1, gcmcparams%nspc
spc = gcmcparams%gcmcspc(j)%spc
fug = eos_getfugacity(gcmcparams%eosparams, spc) ! fug [=] kPa
gcmcparams%gcmcspc(j)%fuglist(i)%fugacity = fug
!** Get the excess chemical potential (B in Adams Notation)
sivolume = volume*1.0e-30*Nav ! convert to m^3/mole
B = Log(fug*1.0e3*sivolume/(Rgas*tk))
!** This is Z/Omega, got from ideal parameters
! wrong: Zig = eos_getConfInteg(gcmcparams%eosparams, j)
! we should be passing the sorbtype , not the index of gcmcmove
Zig = eos_getConfInteg(gcmcparams%eosparams, spc)
ratio = log(Zig)
!** we need (PV)/(RTZ)
B = B - ratio
! Write(0,'(2a,i4,a,f16.2,6f16.10)') __FILE__,": ",__LINE__, &
! " B ", ratio,B,Zig
gcmcparams%gcmcspc(j)%fuglist(i)%B = B
!** Get mu the chemical potential
!** Calculate the DeBroglie wavelength
mass = molecules_getmass(spc)
mass = mass*1.0e-3 ! convert to kg
Lambda = hplanck/Sqrt(twopi*mass/Nav*Rgas/Nav*tk)
murti = (B - Log(sivolume/Nav/Lambda**3))
gcmcparams%gcmcspc(j)%fuglist(i)%murti = murti
! Write(*,*) fug, B, ratio, murti
End Do
End Do
Current code makes pmaps and emaps for only two types of interactions. pmap for LJ, emap for EWALD. There are no combined maps yet. Non-orthogonal cases are still being tested (as of 21 march 2008)
回复此楼