! Program HNCR !-------------------------------------------------------------------- ! ! ! Program to solve the OZ equation with a given closure ! (hnc as currently programmed) for mixtures of charged hard ! sphere fluids, using Allnat's renomalized OZ equation ! (A.R. Allnat, Mol.Phys. 8, 533 (1964). the numerical iteration ! scheme is based on a Newton-Raphson like procedure due ! S. Labik, A. Malijevsky and P. Vonka, Mol.Phys. 56, 709 (1985) ! adapted to fluid mixtures and incorporating the ! iteration strategy for multiple calculations proposed by ! M. Kinoshita and M. Harada, Mol.Phys.,65,599 (1988) following ! the procedure suggested by E.Lomba for simple fluids, Mol.Phys. ! 68, 87 (1989). ! ! Note : ! ! This is a 2-component version. For multiple components ! the program input must be modified and a new subroutine ! for the computation of the initial solution is required. ! ! Note about program units: ! ! This program uses reduced units of length, it is adviceable ! to take the largest ionic diameter as lenght unit. The inverse ! temperature parameter must be then given in the corresponding ! reduced units in accordance with the choice of unit length. ! ! ! Programmed by E. Lomba and J.S. Hoye ! ! For questions related to the program contact E.Lomba ! e-mail: E.Lomba@iqfr.csic.es ! ! Trondheim 19.5.1991 ! ! F90 version: Madrid, 11.08.1999 ! ! ! Libraries required : this program uses subroutines dgetri and dgetrf from ! LAPACK and BLAS computer libraries (or dxml on DEC). ! !---------------------------------------------------------------- ! ! ************************************************************ ! * * ! * input: * ! * * ! * card #1: * ! * * ! * title : character variable describing the comput- * ! * ation * ! * * ! * card #2: * ! * * ! * nr : point number in fourier transform * ! * dr : grid size * ! * tdiff : final error * ! * P pmixt : parameter of mixture in broyles's method * ! * tmax : maximum allowed time for computation (sec)* ! * amx : maximum error to start the change in pmixt* ! * * ! * card #3: * ! * * ! * num : number of points in n-r iterations * ! * (num < or equal numx) * ! * diffnr : error in n-r iterations * ! * * ! * card #4: * ! * * ! * sigma(j) : atomic diameter of species j=1,nsp * ! * It is adviceable to set one of the * ! * diameters to unity (typically the * ! * largest one, sigma) * ! * * ! * card #5: * ! * * ! * qq(j) : charge in e. units of species * ! * j (j=1,nsp) * ! * * ! * card #6: * ! * * ! * nrho,ntr : number of densities and inverse temp. * ! * * ! * card #7: * ! * * ! * densid(j) : total density (sum of all ionic species) * ! * (corresponds to (N/V)*sigma**3 * ! * (j = 1,nrho) * ! * * ! * card #8: * ! * * ! * tk(j) : inverse temperature (e**2/kt*sigma*eps) * ! * eps is the solvent dielectric constant * ! * (j = 1, ntr) * ! * * ! * card #9: * ! * * ! * compsi : if .t. an MSA initial solution is sup- * ! * plied, otherwise is read from disk * ! * * ! ************************************************************ ! * * ! * input files : * ! * hncrin.dat : contains input data * ! * ishnc.dat : contains initial short range * ! * component of r*s(r) * ! * (required if compsi=.f.) * ! * * ! * * ! * output files : * ! * hncrut.dat : Main output file * ! * thermo.dat : Thermodynamics of the systems * ! * hncrg.dat : Ion-Ion correlation functions * ! * ishnc1.dat : short range component of r*s(r) * ! * * ! ************************************************************ ! * * ! * Machine dependent routines : second() CPU time in secs. * ! * See function second in * ! * this code for various options * ! * * ! ************************************************************ ! ! Module section for intercommunication of routines ! Module Params ! ! Specific dimensions are set here ! Integer, Parameter :: NPMX=16384,NUMX=400,NSP=2 Integer, Parameter :: NIT=(NSP*NSP+NSP)/2 Integer, Parameter :: NMT=NIT*NPMX,ND=NIT*NUMX,NMT1=NIT*(NPMX+1) End Module Params Module RQst Use Params ! this module transfers variables to pyhs and msa ! Real (kind=8), Dimension(NPMX) :: R , Q Real (kind=8) :: DQ , DR , CQ1 , CR1 End Module RQst Module FFT ! this module transfers variables to fftsin and fftcos ! and pyhs Use Params Real (kind=8), Dimension(NPMX) :: SINe , COSe , FR , FI Integer :: N , NR , NU End Module FFT ! ! Main program ! Program HNCR Use Params Use RQst Use FFT ! Implicit None Real (kind=8) :: a1 , a1kt , a1sum , a2kt , aa1 , akcor , aklr , & & alfa , alr , alrsum , amx , ankt , asr , chemp , & & chm , chss , chsum , clr , cq Real (kind=8) :: cr , csv , cv , d , dcf , dcfc , & & dcflow , delg , delnr , dels , densid , dethk , & & diff , diffm , diffnr , dnt , dnt2 Real (kind=8) :: dri , eik , ekt , elr , esum , fcnr , fklr , & & fklrii , gamm , gammi , h , hc , hij , & & hij2 , hlow , pi , pmixt Real (kind=8) :: pr , pres , psum , qi , qq , rho , rhoi , & & rhot , ri , rin , rj , rja , rjai , rkd , rkdx2 , & & rkx , rms1 , rms2 Real (kind=8) :: rmusum , rr , s , sclr , sek , sigma , sigx , & & sij , sint , sm1 , sm11 , sm3 , sncj , & & snew , ssm1 , ssum , sum , tc2 , tcc Real (kind=8) :: tcn , tdcf , tdcf0 , tdiff , tfin , th , thij , & & thlr , tim , timer , tk , tmax , tmp , tr , & & trachk , tri , ts , tsn , tsr , tt Real (kind=8) :: vk , vki , work , xc , xcsum , xe , xg , xm , & & xm2 , xp , xp2 , za Integer i , i1 , i2 , ic , ic1 , ic2 , ics , id1 , id2 , idd , & & idn , ik , imin , in1 , in2 , inr , ipk , ipvt , is , it1 Integer it2 , iter , itr , j , js , k , kcut , kwit , l , l1 , & & l2 , mx , ncj , ncmin , ncore , newton , ni Integer npos , nps , nrho , nrmax , & & ntr , ntt , num , nx ! ! ************************************************************ ! * * ! * Most important program parameters : * ! * (to be set in module Params) * ! * npmx = maximum no. of integration points * ! * numx = maximum no. of points in NR iterations * ! * nsp = number of ionic species * ! * nit = number of different interactions * ! * * ! * Note: npmx must be the same in hncr,pyhs,msa,fftcos * ! * and fftsin * ! ************************************************************ ! ! ! ! ! ! Correlation functions (and related quantities) are stored in ! symmetric mode, i.e. ! ! A(I,J) = A(J,I) ---> A((J*J-J)/2+I) (I < or = J) ! Dimension dcf(NPMX,NIT) , tdcf(NPMX,NIT) , ts(NPMX,NIT) , & & h(NPMX,NIT) , tsn(NPMX,NIT) , s(NPMX,NIT) , & & snew(NPMX,NIT) , clr(NPMX,NIT) , sclr(NPMX,NIT) , & & vk(NPMX,NIT) , vki(NPMX,NIT) , th(NPMX,NIT) , & & sint(0:NPMX,NIT) Dimension d(0:NPMX,NIT) , rja(ND,ND) , rjai(ND,ND) , dels(ND) , & & delg(ND) , ipvt(ND) , work(ND) , sigma(NIT) , ncore(NIT)& & , gamm(NIT) , imin(NIT) , cr(NIT) , cq(NIT) , rho(NIT) , & & tcn(NIT) , xg(NIT) , tc2(NIT) , tsr(NIT) , ipg(NIT) , & & ic1(NIT,NIT) , ic2(NIT,NIT) , pr(NIT) , ekt(NIT) , & & rhoi(NIT) , tdcf0(NIT) , diff(NIT) , qq(NSP) , & & gammi(NIT) , tmp(NIT) , chm(NIT) , densid(20) , tk(20) , & & a1kt(NIT) , aklr(NIT) , chemp(NSP) , npos(NSP,NSP) ! ! Save some storage (not too much !) ! Real SECOND Logical ipg , compsi , compj , rgrid Character*80 title Character*1 tit Data pi/3.1415926535898/ Data s , dcf , sint/NMT*0.0 , NMT*0.0 , NMT1*0.0/ Data ipg/NIT*.False./ Data kwit , ntt , compj , rgrid/0 , 0 , .True. , .False./ ! !******************************************************************* ! * ! Program Initialization * ! * !******************************************************************* ! ! ! Initialize time counter ! tt = SECOND() ! Open (3,FILE='hncrut.dat',STATUS='unknown') Open (8,FILE='thermo.dat',STATUS='unknown') Open (9,FILE='hncrin.dat',STATUS='old',ERR=100) ! ! read input data ! Read (9,'(a80)') title Read (9,*) NR , DR , tdiff , pmixt , tmax , amx Read (9,*) num , diffnr Read (9,*) (sigma((j*j+j)/2),j=1,NSP) Read (9,*) (qq(j),j=1,NSP) Read (9,*) nrho , ntr Read (9,*) (densid(j),j=1,nrho) Read (9,*) (tk(j),j=1,ntr) Read (9,*) compsi ! ! check matrix dimensions ! If ( ND.Lt.NIT*num ) Then Write (*,99007) ND/NIT Write (3,99007) ND/NIT Stop Endif ! ! initialize output ! Write (3,99003) title , (j,qq(j),j,sigma((j*j+j)/2),j=1,NSP) Write (8,99003) title , (j,qq(j),j,sigma((j*j+j)/2),j=1,NSP) Write (3,99004) tdiff , pmixt , NR , DR , num , diffnr Write (8,99004) tdiff , pmixt , NR , DR , num , diffnr Write (*,99003) title , (j,qq(j),j,sigma((j*j+j)/2),j=1,NSP) Write (*,99004) tdiff , pmixt , NR , DR , num , diffnr Write (8,99013) (j,j=1,NSP) Write (8,99014) ! ! initialize parameters for ft's ! NU = Int(DLOG(Dble(NR))/DLOG(2.0D0)+0.540) NR = 2**NU N = NR - 1 rin = 1.D0/Dble(4*NR) a1 = pi/Dble(NR) Do i = 1 , NR SINe(i) = DSIN(Dble(i)*a1) COSe(i) = DCOS(Dble(i)*a1) Enddo DQ = pi/(Dble(NR)*DR) CR1 = pi*DR CQ1 = DQ/(8.0D0*pi**2) ! ! initialize control matrices for NR iterations ! mx = 1 Do i = 1 , NSP ipg((i*i+i)/2) = .True. Do j = 1 , i npos(i,j) = mx npos(j,i) = mx If ( mx.Ne.(i*i+i)/2 ) sigma(mx) & & = (sigma((i*i+i)/2)+sigma((j*j+j)/2))/2.D0 nx = 1 Do k = 1 , NSP Do l = 1 , k If ( k.Gt.i ) Then in2 = (k*k-k)/2 + i Else in2 = (i*i-i)/2 + k Endif If ( l.Gt.j ) Then in1 = (l*l-l)/2 + j Else in1 = (j*j-j)/2 + l Endif ic1(mx,nx) = in1 ic2(mx,nx) = in2 nx = nx + 1 Enddo Enddo mx = mx + 1 Enddo Enddo ! ! initialize cutoff radii ! Do j = 1 , NIT ncore(j) = Nint(sigma(j)/DR) ! ! check for adequate grid size ! If ( DABS(ncore(j)*DR-sigma(j)).Gt.1.D-5 ) rgrid = .True. Enddo ! ! If ( rgrid ) Then Write (3,99008) Write (*,99008) Stop Endif Do i = 1 , NR R(i) = Dble(i)*DR Q(i) = Dble(i)*DQ Enddo ! !******************************************************************** ! * ! STEP 1: Enter or compute initial estimate of s(r) * ! * !******************************************************************** ! If ( compsi ) Then ! ! msa initial s(r) (reformulate rho(i) for nsp > 2) ! rhot = densid(1) rho(2) = rhot/(1-qq(2)/qq(1)) rho(1) = -rho(2)*qq(2)/qq(1) Call MSA(s,NR,rho,qq,1/tk(1),sigma(1)) Else ! ! read initial estimate from file ! Open (10,FILE='ishnc.dat',STATUS='old',ERR=200) Read (10,*) rhot , dri , ni , tri , (gammi(j),j=1,NIT) Write (3,'(/'' ** Reading initial solution from disk ...'',/)') Write (*,'(/'' ** Reading initial solution from disk ...'',/)') If ( DABS(DR-dri).Gt.1.E-6 ) Then Write (3,& &'(/'' ** Interpolating S(I,J) to'', '' match input dr **''/)') Write (*,& &'(/'' ** Interpolating S(I,J) to'', '' match input dr **''/)') nrmax = Min(ni,NR-1) Do i = 1 , nrmax Read (10,*) rr , (sint(i,j),j=1,NIT) Enddo nrmax = ni If ( R(NR-1).Lt.rr ) nrmax = NR - 1 Do i = 1 , nrmax inr = Int(R(i)/dri) Do j = 1 , NIT s(i,j) = sint(inr,j) + (sint(inr+1,j)-sint(inr,j))& & *(Dble(i)*DR/dri-inr) Enddo Enddo Else Do i = 1 , Min(ni,NR) Read (10,*) rr , (s(i,j),j=1,NIT) Enddo Endif Endif !********************************* ! loop over densities * !********************************* Do idn = 1 , nrho rhot = densid(idn) ! ! compute ion densities from total density using charge neutrality ! condition (only for nsp = 2) ! rho(2) = rhot/(1-qq(2)/qq(1)) rho(1) = -rho(2)*qq(2)/qq(1) Write (*,'(/ " Calulating for :"/ " rho(1) =",f10.6& &/" rho(2) =",f10.6/" eta =",f10.6/)') rho(1) , rho(2) , & &pi*(rho(1)*sigma(1)**3+rho(2)*sigma(3)**3)/6 rkx = 0.0 Do idd = 1 , NSP rkx = rkx + rho(idd)*qq(idd)**2 Enddo ! !********************************* ! loop over temperatures * !********************************* ! Do itr = 1 , ntr tr = 1/tk(itr) ! ! compute debye inverse screening length ! rkd = DSQRT(4*pi*rkx/tr) ! ! compute adequate yukawa parameter to supress singularity at r=0 ! in short range functions ! If ( rkd.Gt.1 ) Then xe = 4*rkd Else xe = 2.0 Endif xp2 = (xe**2-DSQRT(xe**4-4*(xe*rkd)**2))/2.0 xm2 = (xe**2+DSQRT(xe**4-4*(xe*rkd)**2))/2.0 xp = DSQRT(xp2) xm = DSQRT(xm2) Write (3,99012) 1/tr , rhot Write (*,99012) 1/tr , rhot Write (3,99001) Write (*,99001) mx = 1 ! ! Set density and temperature dependent parameters ! Do id1 = 1 , NSP Do id2 = 1 , id1 ! ! cr(i),cq(i) are density scaled factors for forward and backward ! FFT respectively ! cr(mx) = CR1*DSQRT(rho(id1)*rho(id2)) cq(mx) = CQ1/DSQRT(rho(id1)*rho(id2)) dnt = cr(mx)/CR1 gamm(mx) = -qq(id1)*qq(id2)/tr cv = 4*pi*dnt*gamm(mx)*xe**2 csv = gamm(mx)*xe**2/(xm2-xp2) rkdx2 = (rkd*xe)**2 Do i = 1 , NR qi = Q(i) ri = R(i) ! !****************************************************************** ! * ! STEP 2: Compute and store long range corrections, * ! v(k) and [v(k)]**(-1) elements of the OZ eq. * ! * !****************************************************************** ! ! ! compute the lr behaviour of c(12) and h(12) (note that ! a yukawa term is added so as to make c(12) - clr(12) finite ! at r=0). ! clr(i,mx) = gamm(mx)*(1-DEXP(-xe*ri)) sclr(i,mx) = csv*(DEXP(-xp*ri)-DEXP(-xm*ri)) ! ! compute the chain sum hypervertex v(k) matrix and its inverse ! vk(i,mx) = cv/(qi**4+(qi*xe)**2+rkdx2) vki(i,mx) = -cv/(qi**4+(qi*xe)**2) If ( ipg(mx) ) Then vk(i,mx) = 1 + vk(i,mx) vki(i,mx) = 1 + vki(i,mx) Endif Enddo mx = mx + 1 Enddo Enddo ncmin = NR Do j = 1 , NIT ncmin = Min(ncmin,ncore(j)) ! ! transform series function ! Call FFTSIN(s(1,j),cr(j),ts(1,j)) Enddo ! ! ! begin iterative solution ! iter = 1 ! ! begin broyles cycle ! ! ! ! initiate loop over interactions ! 20 Do j = 1 , NIT ! !****************************************************************** ! * ! STEP 3: Closure relation (c(r) and h(r) computed here) * ! * !****************************************************************** ! ! ! calculate the direct correlation function (hnc closure) ! (lr are subtracted) ! Do i = 1 , ncore(j) - 1 dcf(i,j) = -R(i) - sclr(i,j) - s(i,j) Enddo Do i = ncore(j) , N dcf(i,j) = R(i)*(DEXP((s(i,j)+gamm(j)-clr(i,j)+sclr(i,j))& & /R(i))-1) - sclr(i,j) - s(i,j) Enddo ! ! compute h(12) (including the lr behaviour) ! Do i = 1 , N h(i,j) = (s(i,j)+dcf(i,j)+sclr(i,j))/R(i) Enddo ! !******************************************************************* ! * ! STEP 4: Transform direct correlation function c(12) * ! (the core discontinuity must be treated adequately) * ! * !******************************************************************* ! dcflow = -R(ncore(j)) - sclr(ncore(j),j) - s(ncore(j),j) dcfc = dcf(ncore(j),j) dcf(ncore(j),j) = (dcf(ncore(j),j)+dcflow)/2.0D0 Call FFTSIN(dcf(1,j),cr(j),tdcf(1,j)) dcf(ncore(j),j) = dcfc ! !******************************************************************* ! * ! STEP 5: Initialize functions for NR loop * ! * !******************************************************************* ! ! ! compute d(l) functions (hnc closure) ! hc = h(ncore(j),j) hlow = -1.0D0 - sclr(ncore(j),j) h(ncore(j),j) = (h(ncore(j),j)+hlow)/2.0D0 Call FFTCOS(0.0D0,h(1,j),d(1,j)) sum = 0.0D0 Do i = 1 , N sum = sum + h(i,j) Enddo d(0,j) = 4.D0*sum h(ncore(j),j) = hc ! ! Initial NR guess ! Do i = 1 , num tsn(i,j) = ts(i,j) Enddo Enddo ! ! loop over interactions closed above ! ! !******************************************************************* ! * ! STEP 6: Newton Raphson procedure over j=1,..NUM points * ! * !******************************************************************* ! newton = 0 Do While ( .True. ) ! ! start newton raphson iterations ! ! ! 6.1 *** Check for excess of CPU time ! If ( kwit.Eq.1 ) Then Write (3,'(/" *** convergence not achieved after ", i3," n-r iter and ",i3," broyles iter"/)') newton , iter Write (*,'(/" *** convergence not achieved after ", i3," n-r iter and ",i3," broyles iter"/)') newton , iter Stop Endif ! ! 6.2 *** Compute the jacobian matrix ! Do i = 1 , num Do j = 1 , NIT sum = 0.0D0 Do k = 1 , num ik = IABS(i-k) sum = sum + (d(ik,j)-d(i+k,j))*(tsn(k,j)-ts(k,j)) Enddo tcc = tdcf(i,j) + sum*rin xg(j) = Q(i)*vki(i,j) - tcc tsr(j) = tcc + Q(i)*vk(i,j) Enddo Call SSINV(xg,NSP) Do j = 1 , NIT tsr(j) = xg(j)*Q(i)**2 - tsr(j) dels(num*(j-1)+i) = tsr(j) - tsn(i,j) Enddo ! ! perform the explicit computation and inversion only if ! this is the first system ! If ( compj ) Then Do l1 = 1 , NIT Do l2 = l1 , NIT in1 = ic1(l1,l2) in2 = ic2(l1,l2) rj = -Q(i)**2*xg(in1)*xg(in2) If ( l1.Eq.l2 ) rj = rj + 1.D0 it1 = num*(l1-1) i1 = it1 + i it2 = num*(l2-1) Do k = 1 , num i2 = it2 + k ik = IABS(i-k) rja(i1,i2) = rj*(d(ik,l2)-d(i+k,l2))*rin If ( l1.Ne.l2 ) rja(i+it2,k+it1)& & = rj*(d(ik,l1)-d(i+k,l1))*rin If ( i1.Eq.i2 ) rja(i1,i1) = rja(i1,i1)& & + 1.D0 Enddo Enddo Enddo Endif Enddo ! ! invert jacobian (only if this is the first system) ! If ( compj ) Call INVMAT(rja,ND,NIT*num,rjai,work,ipvt) ! ! 6.2 *** Compute new s(k) using n-r strategy ! ssum = 0.0D0 Do i = 1 , num*NIT sum = 0.0D0 Do j = 1 , num*NIT sum = sum + rjai(i,j)*dels(j) Enddo delg(i) = sum ssum = ssum + sum*sum Enddo delnr = DSQRT(ssum) newton = newton + 1 ntt = ntt + 1 ! ! 6.3 *** Check convergence and time ! timer = SECOND() - tt If ( timer.Gt.tmax ) kwit = 1 If ( delnr.Gt.400 .And. ntt.Gt.1000 ) kwit = 1 ! If ( delnr.Gt.diffnr ) Then Do j = 1 , NIT Do i = 1 , num tsn(i,j) = tsn(i,j) + delg(i+num*(j-1)) Enddo Enddo ! ! if not converged go back to step 6.1 ! Goto 40 Endif ! !******************************************************************* ! * ! STEP 7: Direct OZ equation solution for i=num+1,N * ! * !******************************************************************* ! aa1 = 0.0D0 ! ! check for maximum significant k value ! Do kcut = N , num + 1 , -1 Do j = 1 , NIT aa1 = Max(aa1,DABS(tdcf(kcut,j))) Enddo If ( aa1.Gt.5.0E-4 ) Goto 60 Enddo Goto 60 40 Enddo 60 Do i = num + 1 , N ! ! Here the matrix OZ is solved in k-space ! Do j = 1 , NIT tcc = tdcf(i,j) xg(j) = Q(i)*vki(i,j) - tcc tsr(j) = tcc + Q(i)*vk(i,j) Enddo Call SSINV(xg,NSP) Do j = 1 , NIT tsn(i,j) = xg(j)*Q(i)**2 - tsr(j) Enddo Enddo ! !******************************************************************* ! * ! STEP 8: Transform back to real space s(r) including all k values * ! * !******************************************************************* ! Do j = 1 , NIT Call FFTSIN(tsn(1,j),cq(j),snew(1,j)) Enddo ! !******************************************************************* ! * ! STEP 9: Check for final convergence * ! * !******************************************************************* ! diffm = 0.0D0 Do j = 1 , NIT diff(j) = 0.0D0 Do i = 1 , N diff(j) = diff(j) + (snew(i,j)-s(i,j))**2 Enddo diff(j) = DSQRT(DR*diff(j)) diffm = Max(diffm,diff(j)) Enddo ! !******************************************************************* ! * ! STEP 10: Broyles mixing iterates step in r- and k-space * ! with self-adjusting mixture parameter * ! * !******************************************************************* ! Do j = 1 , NIT If ( diff(j).Gt.amx ) Then alfa = pmixt Else alfa = 1.0D0 - (1.D0-pmixt)*(diff(j)/amx)**2 Endif Do i = 1 , N s(i,j) = s(i,j) + alfa*(snew(i,j)-s(i,j)) ts(i,j) = ts(i,j) + alfa*(tsn(i,j)-ts(i,j)) Enddo Enddo tim = SECOND() - tt Write (3,99002) iter , newton , diffm , tim Write (*,99002) iter , newton , diffm , tim ! ! Branch according to convergence ! If ( diffm.Lt.tdiff ) Then ! ! end broyles cycle ! ! ! set compj so to avoid recomputation of the jacobian inside ! the NR loop ! compj = .False. ! !******************************************************************* ! * ! Compute the jacobian matrix for the last iteration * ! for each system to guarantee stability * ! * !******************************************************************* ! Do i = 1 , num Do j = 1 , NIT sum = 0.0D0 Do k = 1 , num ik = IABS(i-k) sum = sum + (d(ik,j)-d(i+k,j))*(tsn(k,j)-ts(k,j)) Enddo tcc = tdcf(i,j) + sum*rin xg(j) = Q(i)*vki(i,j) - tcc tsr(j) = tcc + Q(i)*vk(i,j) Enddo Call SSINV(xg,NSP) Do j = 1 , NIT tsr(j) = xg(j)*Q(i)**2 - tsr(j) dels(num*(j-1)+i) = tsr(j) - tsn(i,j) Enddo Do l1 = 1 , NIT Do l2 = l1 , NIT in1 = ic1(l1,l2) in2 = ic2(l1,l2) rj = -Q(i)**2*xg(in1)*xg(in2) If ( l1.Eq.l2 ) rj = rj + 1.D0 it1 = num*(l1-1) i1 = it1 + i it2 = num*(l2-1) Do k = 1 , num i2 = it2 + k ik = IABS(i-k) rja(i1,i2) = rj*(d(ik,l2)-d(i+k,l2))*rin If ( l1.Ne.l2 ) rja(i+it2,k+it1) = rj*(d(ik,l1)-d(i+k,l1))*rin If ( i1.Eq.i2 ) rja(i1,i1) = rja(i1,i1) + 1.D0 Enddo Enddo Enddo Enddo ! ! Invert jacobian and sore it in system matrix (rjai) ! Call INVMAT(rja,ND,NIT*num,rjai,work,ipvt) ! !******************************************************************* ! * ! STEP 11: compute thermodynamics * !******************************************************************* ! xcsum = 0.0 Do j = 1 , NIT ncj = ncore(j) dnt2 = (cr(j)/CR1)**2 ic = 0 If ( ipg(j) ) ic = 1 ! ! restore direct correlation function (subtracting the ! coulomb tail for computation of isothermal compressibility) ! and fourier transform of h(r) ! Do i = 1 , N thlr = (vk(i,j)-ic)*Q(i) th(i,j) = ts(i,j) + tdcf(i,j) + thlr dcf(i,j) = dcf(i,j) + clr(i,j) - gamm(j) Enddo ! ! contributions to energy, a/nkt and chemical potentials ! sncj = s(ncj,j) + sclr(ncj,j) - clr(ncj,j) + gamm(j) rms1 = -R(ncj) - sncj ! ! Backward Simpson's rule for r < or = sigma_ij ! Do i = ncj - 1 , 1 , -1 ics = 2*(Mod(ncj-i,2)+1) rms1 = rms1 + ics*dcf(i,j)*R(i) Enddo esum = -R(ncj)*gamm(j)*h(ncj,j) hij2 = (h(ncj,j)*R(ncj))**2 a1sum = -(h(ncj,j)+1)*dcf(ncj,j)*R(ncj) rmusum = R(ncj)*h(ncj,j)*sncj - 2*R(ncj)*dcf(ncj,j) rms2 = dcf(ncj,j)*R(ncj) alrsum = hij2 ! ! Forward Simpson's rule for r > sigma_ij ! Do i = ncj + 1 , N ri = R(i) hij = h(i,j)*ri sij = hij - dcf(i,j) ics = 2*(Mod(i-ncj,2)+1) esum = esum - ics*gamm(j)*hij hij2 = hij*hij a1sum = a1sum - ics*(hij+ri)*dcf(i,j) alrsum = alrsum + ics*hij2 rmusum = rmusum + ics*(hij*sij-2*ri*dcf(i,j)) rms2 = rms2 + ics*dcf(i,j)*ri Enddo ekt(j) = 2*pi*esum*DR/3 + pi*gamm(j)*sigma(j)**2 a1kt(j) = 2*pi*a1sum*DR/3 chm(j) = 2*pi*(rmusum-rms1)*DR/3 xcsum = (2-ic)*(rms1+rms2)*dnt2 + xcsum aklr(j) = 2*pi*alrsum*DR/3 Enddo a2kt = 0.0 ! ! compute integral in k-space for a/nkt (Simpson's rule) ! Do i = 1 , N qi = Q(i) ics = 2*(Mod(i,2)+1) trachk = 0.0 fklr = 0.0 fklrii = 0.0 Do j = 1 , NIT thij = th(i,j)/qi ic = 0 If ( ipg(j) ) Then ic = 1 trachk = trachk + thij fklrii = fklrii + thij*thij Endif xg(j) = ic + thij fklr = fklr + thij*thij Enddo fklr = -(fklr-fklrii) - 0.5*fklrii Call SDET(xg,NSP,dethk) a2kt = a2kt + ics*(DLOG(dethk)-trachk-fklr)*qi*qi Enddo a2kt = -a2kt*DQ/(3.0D0*4*rhot*pi**2) sek = 0.0 psum = 0.0 ankt = 0.0 akcor = 0.0 ! ! add up contributions from each interaction ! l = 1 Do i = 1 , NSP Do j = 1 , i ipk = j/i ic = 2 - ipk fcnr = Dble(ic) a1kt(l) = a1kt(l)*rho(i)*rho(j)/rhot aklr(l) = fcnr*aklr(l)*rho(i)*rho(j)/rhot pr(l) = rho(i)*rho(j)*(h(ncore(l),l)+1)*sigma(l)**3 sek = sek + ic*ekt(l)*rho(i)*rho(j)/rhot ankt = ankt + a1kt(l)*ic psum = psum + pr(l)*ic akcor = akcor + aklr(l) l = l + 1 Enddo Enddo sm1 = 0.0 sm11 = 0.0 sm3 = 0.0 chss = 0.0 ! !******************************************************************* ! * ! Evaluation of final thermodynamic properties * ! * !******************************************************************* ! ! ! add hard core and lr contributions ! Do is = 1 , NSP ssm1 = 0.0 chsum = 0.0 Do js = 1 , NSP nps = npos(is,js) sigx = sigma(nps) sm1 = sm1 + rho(is)*rho(js)*sigx**3 ssm1 = ssm1 + rho(js)*sigx**3 chsum = rho(js)*(chm(nps)+ekt(nps)) + chsum Enddo sm3 = sm3 + (rho(is)*qq(is)*qq(is))**2 sm11 = sm11 + rho(is)*rho(is)*(sigma(npos(is,is)))**3 eik = DEXP(-rkd*R(N)) ! ! 1. Excess chemical potentials (reference ideal gass) ! chemp(is) = 0.5*eik*(0.5*eik-1)*rkd*qq(is)& & **2/tr + chsum + 2*pi*ssm1/3 chss = chss + chemp(is)*rho(is) Enddo alr = (2*rkd**3/(8*pi*rhot)+pi*sm3/(rkd*rhot*tr**2))& & *DEXP(-2*rkd*R(N)) asr = pi*(2*(sm1-sm11)+2*sm11)/(3*rhot) elr = -rkd**3*DEXP(-rkd*R(N))/(8*pi*rhot) ! ! 2. Reduced inverse isothermal compressibility ! xc = 1 - 4*pi*xcsum*DR/(3*rhot) ! ! 3. Excess internal energy (reference hard-sphere system) ! sek = sek + elr ! ! 4. Excess Helmholtz free energy (reference ideal gas) ! ankt = ankt + a2kt + asr + sek + akcor + alr ! ! 5. Compressibility factor from virial equation ! pres = 1.0D0 + 2*pi*psum/rhot/3.0D0 + sek/3.0D0 ! ! 6. Compressibility factor from A/NkT and chemical potential ! za = 1 + chss/rhot - ankt ! tfin = SECOND() - tt ! !******************************************************************* ! * ! STEP 12: Final print out * ! * !******************************************************************* ! ! ! output thermodynamics ! Write (3,99009) ntt , iter , tfin Write (8,99015) tfin , rhot , 1/tr , pres , sek , ankt , & & za , xc , (chemp(j),j=1,NSP) ! ! print out correlation functions ! ! ! store short ranged s(r) and normal s(r) ! Open (2,FILE='ishnc1.dat',STATUS='unknown') Write (2,*) rhot , DR , N , tr , (gamm(j),j=1,NIT) Do i = 1 , N Write (2,'(f10.5,4x,6f12.5)') R(i) , (s(i,j),j=1,NIT) Enddo Close (2) ! ! print out gab(r) ! Open (4,FILE='hncrg.dat',STATUS='unknown') tit = 'g' Write (4,99005) ((tit,j,i,j=1,i),i=1,NSP) Write (4,99006) Do i = ncmin , N Write (4,'(f9.5,4x,6f11.5)') R(i) , (h(i,j)+1,j=1,NIT) Enddo Close (4) Else ! ! iter = iter + 1 Goto 20 Endif Enddo Enddo ! Stop 100 Write (3,99010) Write (*,99010) Stop 200 Write (3,99011) Write (*,99011) ! !******************************************************************* ! * ! Formats and error messages section * ! * !******************************************************************* ! 99001 Format (/2x,'broyles it. n-r it error time'/1x,60('=')) 99002 Format (3x,i3,10x,i3,6x,f8.5,2x,f8.2) 99003 Format (/2x,' program hncr: '/a80,/1x,60('*')& & /6(' qq(',i1,') =',f9.4,' sigma(',i1,')=',f7.4/:)) 99004 Format (' tdiff=',f6.5,' pmixt=',f4.2,' nr=',i4,' dr=',f6.4,& & ' num =',i2,' diffnr=',f6.5//) 99005 Format (/3x,'r',11x,6(a1,'(',2I1,')',6x:)) 99006 Format (1x,70('-')) 99007 Format (/' *** error: num must be < or = ',i3,' ***'/) 99008 Format (/' *** error: grid size must meet the core conditions '/) 99009 Format (/1x,' * end of computation after ',i4,' n-r iter.',& & ' and ',i3,' broyles iter.'/4x,'cpu time :',f7.2,' sec.'/) 99010 Format (/' *** error: file hncrin.dat not found ***'/) 99011 Format (/' *** error: file ishnc.dat for initial s(r) not found',& & '***'/) 99012 Format (/2x,'computing for beta* = ',f9.4,5x,'rho = ',f7.5/1x,& & 60('-')) 99013 Format (3x,'Thermodynamics:'/2x,17('-')/2x,'time',3x,'rho',2x,& & ' beta* ',2x,'z',5x,'u/nkt',4x,'a/nkt',3x,'z(a)',4x,'xc',& & 3x,3('mu/kt(',i1,')',2x:)) 99014 Format (1x,79('*')) 99015 Format (f6.1,3F7.4,6F8.4) ! End Program HNCR Subroutine MSA(Sout,Nr,Rhoe,Qq,Tr,Sigma) Use RQst ! ! programmed by e. lomba, trondheim 23.02.90 ! Computes the MSA solution for equal sized charged hard spheres ! Return value sout(i,j) contains the function h-c. ! Implicit None Real (kind=8) :: b , ck , CQ , CR , cst , ct , e , pi , & & Qq , rho , Rhoe , rk , rlk , s , Sigma , & & sk Real (kind=8) :: Sout , ss , Tr , ts , v , vk , x , xe Integer i , i1 , i2 , mx , n , Nr Real (kind=8) :: k Dimension ts(NPMX) , s(NPMX) , ss(NPMX) , Sout(NPMX,NIT) , Qq(NSP)& & , Rhoe(NSP) cr = cr1 cq = cq1 pi = DACOS(-1.D0) rk = -1/Tr rho = 0.0 Do i = 1 , NSP rho = rho + Rhoe(i)*Qq(i)**2 Enddo e = rho*pi/6.D0 x = DSQRT(24*e/Tr) If ( x.Gt.1 ) Then xe = 4*x Else xe = 2.0 Endif b = (1+x-DSQRT(1+2*x))/x v = 2*b/Tr n = Nr - 1 ! ! compute c(k) and solve oz eq. for ts(k) ! Do i = 1 , Nr k = Q(i) sk = Sin(k) ck = Cos(k) ct = (v*(-k*ck+sk)+v**2*(-k**2*ck+2*k*sk+2*(ck-1))/(4*rk*k))& & /k**2 - rk*ck/k ct = -4*pi*ct ! ! subtract fourier transform of the long range component of c(k) ! rlk = 4*pi*rk*(1/k-k/(k**2+xe**2)) cst = ct - rlk ! ! solve renormalized oz eq. ! vk = k*rho/(k-rho*rlk) ts(i) = k*cst*(vk/rho)**2/(k-cst*vk) - cst Enddo ! ! transform to real space ! Call FFTSIN(ts,CQ,s) Call PYHS(ss,Nr,rho,1.D0) mx = 1 Do i1 = 1 , NSP Do i2 = 1 , i1 Do i = 1 , n Sout(i,mx) = ss(i) + Qq(i1)*Qq(i2)*s(i) Enddo mx = mx + 1 Enddo Enddo End Subroutine MSA ! ! ! Subroutine PYHS(S,Nr,Rho,Dhs) Use RQst ! ! this soubroutine computes the hard spheres function r*s(r), ! dhs is the hard sphere diameter. rho is the hard sphere density ! the approximation used is the percus-yevick. ! Programmed by E.Lomba Trondheim 20.09.90 ! Implicit None Real (kind=8) :: ck , cpi , CQ , CR , Dhs , et , f1 , & & f2 , f3 , pi , qi , qi2 , qk , qr , qr2 , & & qr3 , qr4 Real (kind=8) :: Rho , rl1 , rl2 , rtc , S , sk , tdcf , ts Integer k , n , Nr Dimension ts(NPMX) Dimension S(*) Data pi/3.1415926535898D0/ cr = cr1 cq = cq1 et = pi*Rho*Dhs**3/6.D0 cpi = -4.D0*pi*Dhs*Dhs rl1 = (1.D0+2.D0*et)**2/(1.D0-et)**4 rl2 = -(1.D0+0.5D0*et)**2/(1.D0-et)**4 n = Nr - 1 Do k = 1 , n qk = Q(k) qi = qk*Dhs qi2 = qi*qi sk = DSIN(qi) ck = DCOS(qi) qr = 1.D0/qi qr2 = qr*qr qr3 = qr2*qr qr4 = qr2*qr2 f1 = (sk-qi*ck)*qr2 f2 = (-2.D0+2.D0*qi*sk-qi2*ck+2.D0*ck)*qr3 f3 = (4.D0*(6.D0+qi*(qi2-6.D0)*sk)+(12.D0*qi2-qi2*qi2-24.D0)& & *ck)*qr4*qr tdcf = cpi*(rl1*f1+6.D0*et*rl2*f2+0.5D0*et*rl1*f3) rtc = Rho*tdcf ts(k) = tdcf*rtc/(qk-rtc) Enddo Call FFTSIN(ts,CQ,S) Return End Subroutine PYHS ! ! ! Subroutine FFTSIN(F,C,Tf) Use FFT Implicit None Real (kind=8) :: C , F , Tf , ti , tr Integer i , IBITR , ii , k , k0 , k1 , k2 , l , & & nr1 Dimension F(NPMX) , Tf(NPMX) ! ! this subroutine returns ! ! nr-1 ! tf(i) = c(i)*4*sumk(f(k)*sin(pi*i*k/nr)) ! 1 ! ! f(0) and f(nr) are both set to 0. ! ! see e.o.brigham, "the fast fourier transform", chapter 10. ! Programmed by F.Lado ! k0 = NR/2 nr1 = NR + 1 FR(1) = 0.0D0 F(NR) = 0.0D0 Do i = 1 , k0 FR(i+1) = F(i+i) FI(i) = F(i+i-1) FR(nr1-i) = -FR(i+1) FI(nr1-i) = -FI(i) Enddo ! k = 0 Do l = 1 , NU Do i = 1 , k0 k1 = k + 1 k2 = k1 + k0 tr = FR(k2) ti = FI(k2) FR(k2) = FR(k1) - tr FI(k2) = FI(k1) - ti FR(k1) = FR(k1) + tr FI(k1) = FI(k1) + ti k = k + 1 Enddo Do While ( .True. ) k = k + k0 If ( k.Lt.NR ) Then Do i = 1 , k0 ii = IBITR(k/k0,NU) ii = ii + ii k1 = k + 1 k2 = k1 + k0 tr = FR(k2)*COSe(ii) + FI(k2)*SINe(ii) ti = FI(k2)*COSe(ii) - FR(k2)*SINe(ii) FR(k2) = FR(k1) - tr FI(k2) = FI(k1) - ti FR(k1) = FR(k1) + tr FI(k1) = FI(k1) + ti k = k + 1 Enddo Else k = 0 k0 = k0/2 Goto 100 Endif Enddo 100 Enddo ! Do k = 1 , NR i = IBITR(k-1,NU) + 1 If ( i.Gt.k ) Then tr = FR(k) ti = FI(k) FR(k) = FR(i) FI(k) = FI(i) FR(i) = tr FI(i) = ti Endif Enddo ! Do i = 1 , N Tf(i) = C*(-(FI(i+1)-FI(nr1-i))+(FR(i+1)-FR(nr1-i))*COSe(i)& & +(FI(i+1)+FI(nr1-i))*SINe(i)) Enddo Return End Subroutine FFTSIN ! ! ! Subroutine FFTCOS(F0,F,Tf) Use FFT Implicit None Real (kind=8) :: F , F0 , Tf , ti , tr Integer i , IBITR , ii , k , k0 , k1 , k2 , l , & & nr1 Dimension F(NPMX) , Tf(NPMX) ! ! this subroutine returns ! ! nr-1 ! tf(i) = 4*(f(0)/2+sumk(f(k)*cos(pi*i*k/nr))) ! 1 ! ! f0 is f(0) and f(nr) is set to 0. ! ! see e.o.brigham, "the fast fourier transform", chapter 10. ! Programmed by F.Lado ! k0 = NR/2 nr1 = NR + 1 FR(1) = F0 F(NR) = 0.0D0 Do i = 1 , k0 FR(i+1) = F(i+i) FI(i) = F(i+i-1) FR(nr1-i) = FR(i+1) FI(nr1-i) = FI(i) Enddo ! k = 0 Do l = 1 , NU Do i = 1 , k0 k1 = k + 1 k2 = k1 + k0 tr = FR(k2) ti = FI(k2) FR(k2) = FR(k1) - tr FI(k2) = FI(k1) - ti FR(k1) = FR(k1) + tr FI(k1) = FI(k1) + ti k = k + 1 Enddo Do While ( .True. ) k = k + k0 If ( k.Lt.NR ) Then Do i = 1 , k0 ii = IBITR(k/k0,NU) ii = ii + ii k1 = k + 1 k2 = k1 + k0 tr = FR(k2)*COSe(ii) + FI(k2)*SINe(ii) ti = FI(k2)*COSe(ii) - FR(k2)*SINe(ii) FR(k2) = FR(k1) - tr FI(k2) = FI(k1) - ti FR(k1) = FR(k1) + tr FI(k1) = FI(k1) + ti k = k + 1 Enddo Else k = 0 k0 = k0/2 Goto 100 Endif Enddo 100 Enddo ! Do k = 1 , NR i = IBITR(k-1,NU) + 1 If ( i.Gt.k ) Then tr = FR(k) ti = FI(k) FR(k) = FR(i) FI(k) = FI(i) FR(i) = tr FI(i) = ti Endif Enddo ! Do i = 1 , N Tf(i) = FR(i+1) + FR(nr1-i) + (FI(i+1)+FI(nr1-i))*COSe(i)& & - (FR(i+1)-FR(nr1-i))*SINe(i) Enddo Return End Subroutine FFTCOS ! ! Function IBITR(J,Nu) Implicit None Integer i , IBITR , J , j1 , j2 , Nu ! ! Bit reversal function. Programmed by F.Lado ! j1 = J IBITR = 0 Do i = 1 , Nu j2 = j1/2 IBITR = IBITR + IBITR + (j1-j2-j2) j1 = j2 Enddo Return End Function IBITR ! Subroutine SDET(A,N,Det) Implicit None Real (kind=8) :: A , a1 , Det , sum Integer i , j , j1 , j2 , jm , jp , k , k1 , N Dimension A(N) ! ! this subroutine returns the determinant of the symmetric nxn ! matrix a. symmetric storage mode is used. ! programmed by f.lado ! If ( A(1).Ge.0.0 ) Then Det = A(1) If ( N.Eq.1 ) Return A(1) = DSQRT(A(1)) k1 = 1 Do k = 2 , N A(k1+1) = A(k1+1)/A(1) k1 = k1 + k Enddo j1 = 1 Do j = 2 , N sum = 0.0D0 jm = j - 1 Do i = 1 , jm sum = sum + A(j1+i)**2 Enddo j2 = j1 + j a1 = A(j2) - sum If ( a1.Lt.0.0 ) Goto 100 Det = Det*a1 A(j2) = DSQRT(a1) If ( j.Ne.N ) Then jp = j + 1 k1 = j2 Do k = jp , N sum = 0.0 Do i = 1 , jm sum = sum + A(j1+i)*A(k1+i) Enddo A(k1+j) = (A(k1+j)-sum)/A(j2) k1 = k1 + k Enddo Endif j1 = j2 Enddo Return Endif 100 Print * , ' *** error in routine sdet ***' End Subroutine SDET ! ! Subroutine SPROD(A,B,C,N) Implicit None Real (kind=8) :: A , B , C , sum Integer i , i1 , is , j , j1 , js , k , ks , N Dimension A(10) , B(10) , C(10) ! ! this subroutine returns the matrix product c = ab for symmetric ! nxn matrices a, b, and c. symmetric storage mode is used. ! ! adapted from the ibm "scientific subroutine package" by F.Lado. ! js = 0 Do j = 1 , N is = 0 Do i = 1 , j sum = 0.0D0 Do k = 1 , i sum = sum + A(is+k)*B(js+k) Enddo ks = is + i If ( i.Ne.j ) Then i1 = i + 1 Do k = i1 , j sum = sum + A(ks+i)*B(js+k) ks = ks + k Enddo Endif If ( j.Ne.N ) Then j1 = j + 1 Do k = j1 , N sum = sum + A(ks+i)*B(ks+j) ks = ks + k Enddo Endif C(js+i) = sum is = is + i Enddo js = js + j Enddo Return End Subroutine SPROD ! ! ! Subroutine SSINV(Ai,N) Implicit None Real (kind=8) :: a , Ai , big , epsiln , p , pivot1 , q , test Integer i , j , k , km1 , kp1 , l , N Integer r Dimension Ai(10) , a(6,6) , p(100) , q(100) , r(100) Data epsiln/1.D-20/ ! ! Inversion of a symmetric matrix. Symmetric storage mode is used. ! Orignal source code from CERN computer library. ! If ( N.Lt.1 ) Then Elseif ( N.Eq.1 ) Then Ai(1) = 1.D0/Ai(1) Return Else l = 1 Do i = 1 , N Do j = 1 , i a(j,i) = Ai(l) l = l + 1 Enddo Enddo ! ! construct truth table ! Do i = 1 , N r(i) = 1 Enddo ! ! begin programme ! Do i = 1 , N k = 0 ! ! search for pivot ! big = 0. Do j = 1 , N test = DABS(a(j,j)) If ( test.Gt.big ) Then If ( r(j).Lt.0 ) Goto 100 If ( r(j).Ne.0 ) Then big = test k = j Endif Endif Enddo ! ! test for zero matrix ! If ( k.Le.0 ) Goto 100 ! ! test for linearity ! If ( i.Eq.1 ) pivot1 = a(k,k) If ( Abs(a(k,k)/pivot1).Lt.epsiln ) Goto 100 ! ! preparation for elimination step1 ! r(k) = 0 q(k) = 1./a(k,k) p(k) = 1. a(k,k) = 0.0 kp1 = k + 1 km1 = k - 1 If ( km1.Lt.0 ) Goto 100 If ( km1.Ne.0 ) Then Do j = 1 , km1 p(j) = a(j,k) q(j) = a(j,k)*q(k) If ( r(j).Lt.0 ) Goto 100 If ( r(j).Ne.0 ) q(j) = -q(j) a(j,k) = 0. Enddo Endif If ( k.Lt.N ) Then Do j = kp1 , N p(j) = a(k,j) q(j) = -a(k,j)*q(k) If ( r(j).Lt.0 ) Goto 100 If ( r(j).Eq.0 ) p(j) = -p(j) a(k,j) = 0.0 Enddo Elseif ( k.Ne.N ) Then Goto 100 Endif ! ! elimination proper ! Do j = 1 , N Do k = j , N a(j,k) = a(j,k) + p(j)*q(k) Enddo Enddo Enddo ! ! elements of left diagonal ! l = 1 Do j = 1 , N Do k = 1 , j Ai(l) = a(k,j) l = l + 1 Enddo Enddo Return Endif ! ! failure return ! 100 Print * , 'failure return from subroutine ssinv' End Subroutine ssinv ! ! ! Subroutine INVMAT(A,Lda,N,B,Work,Ipvt) ! ! Computes the inverse of a matrix ! using routines from LAPACK and BLAS ! Programmed by E.Lomba as inteface to dgetrf and dgtri routines !***description ! ! on entry ! ! a real*8(lda, n) ! the output from sgeco or sgefa. ! ! lda integer ! the leading dimension of the array a . ! ! n integer ! the order of the matrix a . ! ! ipvt integer(n) ! the pivot vector from sgeco or sgefa. ! ! work real(n) ! work vector. contents destroyed. ! ! on return ! ! b inverse of original a matrix ! Implicit None Real (kind=8) :: A , B , det , Work Integer i , info , Ipvt , j , Lda , N Dimension A(Lda,1) , B(Lda,1) , Work(N) , Ipvt(N) , det(2) Do i = 1 , N Do j = 1 , N B(i,j) = A(i,j) Enddo Enddo Call DGETRF(N,N,B,Lda,Ipvt,info) If ( info.Eq.0 ) Then Call DGETRI(N,B,Lda,Ipvt,Work,N,info) Return Endif Write (*,'(/'' *** error from invmat: a is singular ***'')') Stop End Subroutine INVMAT