1C> \ingroup task
2C> @{
3       function task_rism(rtdb)
4       implicit none
5#include "errquit.fh"
6#include "rtdb.fh"
7#include "mafdecls.fh"
8#include "inp.fh"
9#include "util.fh"
10#include "global.fh"
11       integer rtdb
12       logical  task_rism
13       integer n
14       logical rism_util_power_2
15       external rism_util_power_2
16c
17       call rism_print_header()
18c
19c      load data and allocate memory
20c      -----------------------------
21       call rism_prepare(rtdb)
22#ifdef RISM_DEBUG
23       call rism_message("finished with rism_prepare")
24#endif
25       call rism_print_params(rtdb)
26       call rism_print_solute_configuration(rtdb)
27c       stop
28c
29c      main rism routine call
30c      ----------------------
31       call rism_message("calling rism wrapper")
32       call rism_wrapper(rtdb)
33       call rism_cleanup()
34       task_rism = .true.
35       call rism_message("completed task rism")
36       return
37      end
38C> @}
39c
40       subroutine rism_wrapper(rtdb)
41       implicit none
42#include "errquit.fh"
43#include "rtdb.fh"
44#include "mafdecls.fh"
45#include "inp.fh"
46#include "util.fh"
47#include "global.fh"
48      integer rtdb
49      integer nu,nv,nvv,ngr
50      integer icl,icr
51      integer i_rgrid,i_kgrid
52      integer i_tv,i_den,i_isv,i_mv
53      integer i_xv,i_yv,i_zv
54      integer i_ims
55      logical result
56      character*32 sname,dname,pname
57      double precision t,tau,lambd,tol,solperm
58      integer dd
59      integer i_sgvv,i_epsvv,i_qvv
60      integer i_sigu,i_epsiu,i_qqu,i_wu
61      character*72 rdffile,tag
62      logical okspace
63c
64      pname = "rism_wrapper"
65c
66c
67c     solute data
68c     --------------
69      sname = "solute"
70      dname = "natoms"
71      call db_data_get_int(sname,dname,1,nu,result)
72      if(.not.result)
73     >  call errquit(pname//"cannot get"//sname(1:inp_strlen(sname))
74     >               //dname(1:inp_strlen(dname)),0,0)
75
76      dname = "sigma"
77      call db_data_get_index(sname,dname,i_sigu,result)
78      if(.not.result)
79     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
80     >               //dname(1:inp_strlen(dname)),0,0)
81
82      dname = "epsilon"
83      call db_data_get_index(sname,dname,i_epsiu,result)
84      if(.not.result)
85     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
86     >               //dname(1:inp_strlen(dname)),0,0)
87
88      dname = "charge"
89      call db_data_get_index(sname,dname,i_qqu,result)
90      if(.not.result)
91     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
92     >               //dname(1:inp_strlen(dname)),0,0)
93
94      dname = "struct_factor"
95      call db_data_get_index(sname,dname,i_wu,result)
96      if(.not.result)
97     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
98     >               //dname(1:inp_strlen(dname)),0,0)
99
100
101#ifdef RISM_DEBUG
102      call rism_message("finished with solute data")
103#endif
104c
105c     solvent data
106c     --------------
107      sname = "solvent"
108      dname = "natoms"
109      call db_data_get_int(sname,dname,1,nv,result)
110      if(.not.result)
111     >  call errquit(pname//"cannot get"//sname(1:inp_strlen(sname))
112     >               //dname(1:inp_strlen(dname)),0,0)
113
114      dname = "natoms_reduced"
115      call db_data_get_int(sname,dname,1,nvv,result)
116      if(.not.result)
117     >  call errquit(pname//"cannot get"//sname(1:inp_strlen(sname))
118     >               //dname(1:inp_strlen(dname)),0,0)
119
120      dname = "atom_name"
121      call db_data_get_index(sname,dname,i_tv,result)
122      if(.not.result)
123     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
124     >               //dname(1:inp_strlen(dname)),0,0)
125
126      dname = "residue_index"
127      call db_data_get_index(sname,dname,i_isv,result)
128      if(.not.result)
129     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
130     >               //dname(1:inp_strlen(dname)),0,0)
131
132      dname = "density"
133      call db_data_get_index(sname,dname,i_den,result)
134      if(.not.result)
135     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
136     >               //dname(1:inp_strlen(dname)),0,0)
137
138      dname = "multiplicity"
139      call db_data_get_index(sname,dname,i_mv,result)
140      if(.not.result)
141     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
142     >               //dname(1:inp_strlen(dname)),0,0)
143
144      dname = "xcoord"
145      call db_data_get_index(sname,dname,i_xv,result)
146      if(.not.result)
147     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
148     >               //dname(1:inp_strlen(dname)),0,0)
149
150      dname = "ycoord"
151      call db_data_get_index(sname,dname,i_yv,result)
152      if(.not.result)
153     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
154     >               //dname(1:inp_strlen(dname)),0,0)
155
156      dname = "zcoord"
157      call db_data_get_index(sname,dname,i_zv,result)
158      if(.not.result)
159     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
160     >               //dname(1:inp_strlen(dname)),0,0)
161
162      dname = "map_reduced"
163      call db_data_get_index(sname,dname,i_ims,result)
164      if(.not.result)
165     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
166     >               //dname(1:inp_strlen(dname)),0,0)
167
168
169      dname = "sigma_reduced"
170      call db_data_get_index(sname,dname,i_sgvv,result)
171      if(.not.result)
172     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
173     >               //dname(1:inp_strlen(dname)),0,0)
174
175      dname = "epsilon_reduced"
176      call db_data_get_index(sname,dname,i_epsvv,result)
177      if(.not.result)
178     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
179     >               //dname(1:inp_strlen(dname)),0,0)
180
181      dname = "charge_reduced"
182      call db_data_get_index(sname,dname,i_qvv,result)
183      if(.not.result)
184     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
185     >               //dname(1:inp_strlen(dname)),0,0)
186      call rism_message("rism_wrapper 1")
187c
188c     grid data
189c     --------------
190      sname = "grid"
191      dname = "npoints"
192      call db_data_get_int(sname,dname,1,ngr,result)
193      if(.not.result)
194     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
195     >               //dname(1:inp_strlen(dname)),0,0)
196
197      dname = "rgrid"
198      call db_data_get_index(sname,dname,i_rgrid,result)
199      if(.not.result)
200     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
201     >               //dname(1:inp_strlen(dname)),0,0)
202
203      dname = "okspace"
204      call db_data_get_log(sname,dname,1,okspace,result)
205      if(.not.result)
206     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
207     >               //dname(1:inp_strlen(dname)),0,0)
208      dname = "kgrid"
209      call db_data_get_index(sname,dname,i_kgrid,result)
210      if(.not.result)
211     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
212     >               //dname(1:inp_strlen(dname)),0,0)
213      call rism_message("rism_wrapper 2")
214c
215c     parameters
216c     --------------
217      sname = "parameters"
218      dname = "closure"
219      call db_data_get_int(sname,dname,1,icl,result)
220      if(.not.result)
221     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
222     >               //dname(1:inp_strlen(dname)),0,0)
223
224      dname = "vdw_rule"
225      call db_data_get_int(sname,dname,1,icr,result)
226      if(.not.result)
227     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
228     >               //dname(1:inp_strlen(dname)),0,0)
229
230      dname = "solvent_permittivity"
231      call db_data_get_dbl(sname,dname,1,solperm,result)
232      if(.not.result)
233     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
234     >               //dname(1:inp_strlen(dname)),0,0)
235
236      dname = "tau"
237      call db_data_get_dbl(sname,dname,1,tau,result)
238      if(.not.result)
239     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
240     >               //dname(1:inp_strlen(dname)),0,0)
241
242
243      dname = "tolerance"
244      call db_data_get_dbl(sname,dname,1,tol,result)
245      if(.not.result)
246     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
247     >               //dname(1:inp_strlen(dname)),0,0)
248
249      dname = "mixing"
250      call db_data_get_dbl(sname,dname,1,lambd,result)
251      if(.not.result)
252     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
253     >               //dname(1:inp_strlen(dname)),0,0)
254
255      dname = "temperature"
256      call db_data_get_dbl(sname,dname,1,t,result)
257      if(.not.result)
258     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
259     >               //dname(1:inp_strlen(dname)),0,0)
260
261      dname = "diis"
262      call db_data_get_int(sname,dname,1,dd,result)
263      if(.not.result)
264     >  call errquit(pname//"cannot get "//sname(1:inp_strlen(sname))
265     >               //dname(1:inp_strlen(dname)),0,0)
266c
267c     get filename for solvent g(r)
268c     -----------------------------
269      tag = "rism:solvent:rdf"
270      if(.not.rtdb_cget(rtdb,tag,1,rdffile))
271     >  call errquit("cannot get "//tag,0,0)
272      call rism_message("rism_wrapper 3")
273
274       call rism_message("getting ready for main rism")
275
276       call rism(rtdb,rdffile,nu,nv,nvv,ngr,icl,
277     +           dbl_mb(i_kgrid),dbl_mb(i_rgrid),
278     +           byte_mb(i_tv),int_mb(i_isv),dbl_mb(i_den),
279     +           int_mb(i_mv),
280     +           dbl_mb(i_xv),dbl_mb(i_yv),dbl_mb(i_zv),
281     +           icr,int_mb(i_ims),tau,solperm,
282     +           dbl_mb(i_sgvv),dbl_mb(i_epsvv),
283     +           dbl_mb(i_qvv),dbl_mb(i_sigu),
284     +           dbl_mb(i_epsiu),dbl_mb(i_qqu),
285     +           t,lambd,
286     +           tol,dbl_mb(i_wu),dd,okspace)
287
288c       call rism(nu,nv,nvv,ngr,icl,kgrid,rgrid,tv,isv,den,mv,xv,yv,zv,
289c     * icr,ims,tau,solperm,sgvv,epsvv,qvv,sigu,epsiu,qqu,t,lambd,
290c     * tol,wu,dd)
291
292       return
293      end
294
295c creat susceptibility of solvent
296c
297c  creat solute-solvent potentials
298c
299       subroutine potcreat(nvv,nu,ngr,icr,kgrid,rgrid,tau,solperm,
300     * sgvv,epsvv,qvv,sigu,epsiu,qqu,plj,ul)
301       implicit none
302       integer  nu,nv, ngr,nvv
303       real*8 sigu(1:nu),epsiu(1:nu),qqu(1:nu)
304       real*8 sgvv(1:nvv),epsvv(1:nvv),qvv(1:nvv)
305       integer  i, j1,j2,icr
306       real*8 r1, r2, dr, dk, rgrid(1:ngr),kgrid(1:ngr),pi
307       real*8 sigvv(1:nu,1:nvv), epsivv(1:nu,1:nvv)
308       real*8 ul(1:nu,1:nvv,1:ngr),plj(1:nu,1:nvv,1:ngr)
309       real*8 nav,tau, echar,solperm,uthre,ck
310       external util_erf
311       double precision util_erf
312c      constants for potential
313       pi=2*asin(1.0)
314       nav=6.02214179e+23
315       echar=4.803e-10
316       ck=nav*echar**2/solperm*0.01
317       uthre=1000
318c
319c calculation lj parameters
320       call combrule(nu,nvv,icr,sgvv,sigu,epsiu,epsvv,sigvv,epsivv)
321c
322c plj(j1,j2,r) lj potential+short part from coulomb
323c ul(j1,j2,k) is the fourier trasnform of long range interactions
324       do i=1,ngr
325        do j1=1,nu
326         do j2=1,nvv
327          ul(j1,j2,i)=4*pi*ck*qqu(j1)*qvv(j2)*exp(-kgrid(i)**2/4/tau**2)
328          ul(j1,j2,i)=ul(j1,j2,i)/kgrid(i)**2
329          plj(j1,j2,i)=4*epsivv(j1,j2)*
330     *    ((sigvv(j1,j2)/rgrid(i))**(12)-(sigvv(j1,j2)/rgrid(i))**(6))
331          plj(j1,j2,i)=plj(j1,j2,i)+ck
332     *    *qqu(j1)*qvv(j2)/rgrid(i)*(1-util_erf(tau*rgrid(i)))
333          if (plj(j1,j2,i).ge. uthre) then
334           plj(j1,j2,i)=uthre
335          endif
336         enddo
337        enddo
338       enddo
339!      do i=1,40
340!       print*, (plj(1,j1,i), j1=1,4)
341!      enddo
342       return
343       end subroutine
344c
345c   calculations of parameters for solute-solvent interactions
346c
347       subroutine combrule(nu,nvv,icr,sgvv,sigu,
348     * epsiu,epsvv,sigvv,epsivv)
349       implicit none
350       integer  nu, nv, icr,nvv
351       real*8 sigu(1:nu), epsiu(1:nu), sgvv(1:nvv), epsvv(1:nvv)
352       real*8 sigvv(1:nu,1:nvv), epsivv(1:nu,1:nvv)
353       integer  i, j1,j2
354       do j1=1,nu
355        do j2=1,nvv
356         epsivv(j1,j2)=(epsiu(j1)*epsvv(j2))**(1.0/2)
357         if (icr.eq.1) then
358          sigvv(j1,j2)=(sigu(j1)+sgvv(j2))/2
359          else
360          sigvv(j1,j2)=(sigu(j1)*sgvv(j2))**(1.0/2)
361         endif
362!        print *, sigvv(j1,j2), epsivv(j1,j2), j1,j2
363        enddo
364       enddo
365       return
366       end subroutine
367c
368c site-site ornstein-zernike in k-space
369c
370       subroutine ssoz(nvv,nu,ngr,nv,ims,mv,wu,c3,chi,hu)
371       implicit none
372       integer  nu, nvv, nv, ngr, i,j1,j,k1,ims(1:nv),mv(1:nv)
373       real*8 rgrid(1:ngr),kgrid(1:ngr)
374       real*8 pi,ct
375       real*8 c3(1:nu,1:nvv,1:ngr),hu(1:nu,1:nvv,1:ngr)
376       real*8 ctt(1:nu,1:nvv,1:ngr),h1(1:nu,1:nvv,1:ngr)
377       real*8 wu(1:nu,1:nu,1:ngr), chi(1:nvv,1:nvv,1:ngr)
378c      calculations of right product of matrix via summation over solvent sites
379       do i=1,ngr
380        do j1=1,nu
381         do j=1,nvv
382          ct=0
383          do k1=1,nvv
384           ct=c3(j1,k1,i)*chi(k1,j,i)+ct
385          enddo
386          ctt(j1,j,i)=ct
387         enddo
388        enddo
389       enddo
390!      do i=1,ngr,40
391!       print*, (ctt(1,j,i), j=1,4)
392!      enddo
393c      calculations of left product of matrix via summation over solute sites
394       do i=1,ngr
395        do j1=1,nu
396         do j=1,nvv
397          ct=0
398          do k1=1,nu
399           ct=wu(j1,k1,i)*ctt(k1,j,i)+ct
400          enddo
401          hu(j1,j,i)=ct
402          h1(j1,j,i)=hu(j1,j,i)
403         enddo
404        enddo
405       enddo
406       do i=1,ngr
407        do j1=1,nu
408         do j=1,nv
409          hu(j1,ims(j),i)=h1(j1,ims(j),i)/mv(j)
410         enddo
411        enddo
412       enddo
413!      do i=1,ngr,40
414!       print*, (hu(1,j,i),h1(1,j,i), j=1,4)
415!      enddo
416       return
417       end subroutine
418c
419c  subroutine for closure
420c
421       subroutine clos(nvv,nu,ngr,t,plj,cr,gam,icl)
422       implicit none
423       integer  nu, nvv, ngr, i,j1,j,icl
424       real*8 plj(1:nu,1:nvv,1:ngr)
425       real*8 t, kb,pi,lexp
426       real*8 cr(1:nu,1:nvv,1:ngr),gam(1:nu,1:nvv,1:ngr)
427       kb=8.13441e-3
428       do i=1,ngr
429        do j=1,nu
430         do j1=1,nvv
431          if(icl.eq.1) then
432           cr(j,j1,i)=exp(-1.0/kb/t*plj(j,j1,i)+gam(j,j1,i))
433     *     -gam(j,j1,i)-1
434          endif
435          if(icl.eq.2) then
436           cr(j,j1,i)=lexp(-1.0/kb/t*plj(j,j1,i)+gam(j,j1,i))
437     *     -gam(j,j1,i)-1
438          endif
439         enddo
440        enddo
441       enddo
442       return
443       end subroutine
444c
445       real*8 function lexp(x)
446       real*8 x
447       if(x.le.0.0)then
448        lexp=exp(x)
449        else
450        lexp=1+x
451       endif
452       return
453       end
454c
455ccc
456c
457       subroutine rism(rtdb,rdffile,nu,nv,nvv,ngr,
458     * icl,kgrid,rgrid,tv,isv,
459     * den,mv,xv,yv,zv,
460     * icr,ims,tau,solperm,sgvv,epsvv,qvv,sigu,epsiu,qqu,t,
461     * lambd,tol,wu,dd,okspace)
462       implicit none
463#include "rtdb.fh"
464#include "mafdecls.fh"
465c      parameters
466       integer rtdb
467       integer kd,dd
468       character*(*) rdffile
469       integer  nu, nv,nd,nvv, ngr, i,j1,j,k1,icl,k
470       real*8 rgrid(1:ngr),kgrid(1:ngr)
471       real*8 pi,dk,t,tau,solperm,tol,lambd,del0,kb,dr
472c      solute
473       real*8 wu(1:nu,1:nu,1:ngr), chi(1:nvv,1:nvv,1:ngr)
474       real*8 sigu(1:nu),epsiu(1:nu),qqu(1:nu)
475c      solvent
476       integer  icr,isv(1:nv),mv(1:nv),ims(1:nv)
477       real*8 den(1:nv),xv(1:nv),yv(1:nv),zv(1:nv)
478       real*8 wv(1:nvv,1:nvv,1:ngr)
479       real*8 sgvv(1:nvv),epsvv(1:nvv),qvv(1:nvv)
480       character (4) tv(1:nv)
481       logical okspace
482c      potential
483       real*8 plj(1:nu,1:nvv,1:ngr), ul(1:nu,1:nvv,1:ngr)
484c      functions
485       real*8 cr(1:nu,1:nvv,1:ngr),c2(1:ngr),tt(1:nu,1:nvv)
486       real*8 ht(1:ngr),c3(1:nu,1:nvv,1:ngr),cf3(1:nu,1:nvv,1:ngr)
487       real*8 cold(1:nu,1:nvv,1:ngr),cnew(1:nu,1:nvv,1:ngr)
488       real*8 gfold(1:nu,1:nvv,1:ngr),hu(1:nu,1:nvv,1:ngr)
489       real*8 gold(1:nu,1:nvv,1:ngr), gnew(1:nu,1:nvv,1:ngr)
490       real*8 hold(1:nu,1:nvv,1:ngr), hnew(1:nu,1:nvv,1:ngr)
491       real*8 del(1:nu,1:nvv),mmaxi,mmax(1:nu),mmmax
492c      diis functions and counts
493       real*8 ggo(1:dd,1:nu,1:nvv,1:ngr),dgg(1:dd,1:nu,1:nvv,1:ngr)
494       double precision muh,mugf
495       write(76,*) nu,nv,nvv,ngr,icl,icr,tau,solperm
496       write(76,*) T,lambd,tol
497
498c
499       pi=2*asin(1.0)
500c      nd total number of different rdfs
501       nd=(nvv+1)*nvv/2
502c      chi(i,j,k) solvent susceptibility
503c
504c      WORK HERE
505       call chicreat(rdffile,nd,nv,ngr,ims,nvv,
506     +               kgrid,tv,isv,den,mv,xv,yv,zv,chi,
507     +               qvv,tau,okspace)
508c
509c       call chicreat(rdffile,nd,nv,ngr,ims,nvv,
510c     +               kgrid,tv,isv,den,mv,xv,yv,zv,chi)
511
512c
513c      plj(j1,j2,r) lj potential+short part from coulomb
514c      ul(j1,j2,k) is the fourier trasnform of long range interactions
515       call potcreat(nvv,nu,ngr,icr,kgrid,rgrid,tau,solperm,
516     * sgvv,epsvv,qvv,sigu,epsiu,qqu,plj,ul)
517c
518cc     initial value of c_short(r)
519       kb=8.3144e-3
520       do i=1,ngr
521        do j=1,nu
522         do j1=1,nvv
523          gold(j,j1,i)=0
524          cold(j,j1,i)=0
525         enddo
526        enddo
527       enddo
528c      starting counts for diis (kd) and for iterations k1
529       k1=1
530       kd=1
531       write(*,*) "starting iterations"
5321      continue
533c      calculating c_new(r) and h_new(r) from closure
534       call clos(nvv,nu,ngr,t,plj,cnew,gold,icl)
535       do j=1,nu
536        do j1=1,nvv
537         do i=1,ngr
538          hold(j,j1,i)=gold(j,j1,i)+cnew(j,j1,i)
539         enddo
540        enddo
541       enddo
542c      fast fourier transform
543       dr=rgrid(2)-rgrid(1)
544       do j=1,nu
545        do j1=1,nvv
546         c2(1)=0
547         do i=1,ngr-1
548          c2(i+1)=cnew(j,j1,i)*rgrid(i)
549         enddo
550         call sinft(c2,ngr)
551c        normalization of sin-fft with excluding the zeropoint (x=0)
552         do i=1,ngr-1
553          cf3(j,j1,i)=4*pi*c2(i+1)/kgrid(i)*dr
554         enddo
555         cf3(j,j1,ngr)=cf3(j,j1,ngr-1)
556        enddo
557       enddo
558       dk=kgrid(2)-kgrid(1)
559       pi=2*asin(1.0)
560c       adding the long-range potential to c_short
561       do i=1,ngr
562        do j=1,nu
563         do j1=1,nvv
564          cf3(j,j1,i)=cf3(j,j1,i)-1.0/kb/t*ul(j,j1,i)
565         enddo
566        enddo
567       enddo
568c      calculations of fourier transform of h(i,j,k) by
569c      site-site ornstein-zerinike equations
570       call ssoz(nvv,nu,ngr,nv,ims,mv,wu,cf3,chi,hu)
571c      inverse fourier transform of h(i,j,k) & calculations of gamma(i,i,r)
572       do j=1,nu
573        do j1=1,nvv
574         ht(1)=0
575         do i=1,ngr-1
576          ht(i+1)=hu(j,j1,i)*kgrid(i)
577         enddo
578         call sinft(ht,ngr)
579         do i=1,ngr-1
580          hnew(j,j1,i)=ht(i+1)*dk/2/rgrid(i)/pi**2
581          gnew(j,j1,i)=hnew(j,j1,i)-cnew(j,j1,i)
582         enddo
583         gnew(j,j1,ngr)=gnew(j,j1,ngr-1)
584         hnew(j,j1,ngr)=hnew(j,j1,ngr-1)
585        enddo
586       enddo
587c      intialization del
588       do j=1,nu
589        do j1=1,nvv
590         del(j,j1)=0
591        enddo
592       mmax(j)=0
593       enddo
594       mmaxi=0
595c      evaluation of accuracy
596       do j=1,nu
597        do j1=1,nvv
598         do i=1,ngr
599          del(j,j1)=del(j,j1)+(gnew(j,j1,i)-gold(j,j1,i))**2
600         enddo
601         mmax(j)=del(j,j1)+mmax(j)
602        enddo
603        mmaxi=mmaxi+mmax(j)
604       enddo
605c      normalization
606       del0=(mmaxi/ngr/nu/nvv)**(1.0/2)
607c      diis procedure
608       call  diis(nu,nvv,ngr,lambd,kd,dd,gold,gnew,ggo,dgg)
609c      old functions as new functions for new cycle
610       do j=1,nu
611        do j1=1,nvv
612         do i=1,ngr
613          gold(j,j1,i)=gnew(j,j1,i)
614          cold(j,j1,i)=cnew(j,j1,i)
615         enddo
616        enddo
617       enddo
618       k1=k1+1
619       write(*,*) k1,kd,del0
620       if(k1.ge.500) goto 2
621        if(del0.ge.tol) goto 1
622 2     continue
623       call thermo(nu,nv,nvv,ngr,icl,ims,den,t,rgrid,qqu,qvv,
624     * hnew,cnew,solperm,tau,muh,mugf)
625       if (.not. rtdb_put(rtdb, 'rism:energy', mt_dbl, 1, mugf))
626     $           call errquit('failed rism',0,0)
627       if (.not. rtdb_put(rtdb, 'rism:muh', mt_dbl, 1, muh))
628     $           call errquit('failed rism',0,0)
629       if (.not. rtdb_put(rtdb, 'rism:mugf', mt_dbl, 1, mugf))
630     $           call errquit('failed rism',0,0)
631
632
633       return
634       end subroutine
635cc
636c
637c      subroutine for evaluation diis residues dgg and matrix coefficients ss
638c
639       subroutine diis(nu,nvv,ngr,lam,k1,dd,hold,hnew,ggo,dgg)
640       implicit none
641       integer k1,dd,jd,jd1
642       real*8 ggo(1:dd,1:nu,1:nvv,1:ngr),dgg(1:dd,1:nu,1:nvv,1:ngr)
643       real*8 ss(1:dd+1,1:dd+1),hold(1:nu,1:nvv,1:ngr)
644       real*8 s0(1:dd+1),s1,lam,hnew(1:nu,1:nvv,1:ngr)
645       integer j,j1,i,nu,nvv,ngr
646       integer ipiv(1:dd+1),info
647c      initialization of overlap matrix for diis
648       do jd=1,dd
649        do jd1=1,dd
650         ss(jd,jd1)=0
651        enddo
652        s0(jd)=0
653        ss(jd,dd+1)=-1
654        ss(dd+1,jd)=-1
655       enddo
656       ss(dd+1,dd+1)=0
657       s0(dd+1)=-1
658c      calculation of basis
659       do j=1,nu
660        do j1=1,nvv
661         do i=1,ngr
662          ggo(k1,j,j1,i)=hnew(j,j1,i)
663          dgg(k1,j,j1,i)=hnew(j,j1,i)-hold(j,j1,i)
664c         calculation of elements of overlaping  matrix
665          if(k1.eq.dd) then
666           do jd=1,dd
667            do jd1=jd,dd
668             ss(jd,jd1)=(dgg(jd,j,j1,i)*dgg(jd1,j,j1,i))/ngr+ss(jd,jd1)
669c            symmerization
670             ss(jd1,jd)=ss(jd,jd1)
671            enddo
672           enddo
673          endif
674         enddo
675        enddo
676       enddo
677c      solution of diis equation and change of baisis
678       if(k1.eq.dd) then
679        do jd=1,dd
680        enddo
681        call dgesv(dd+1,1,ss,dd+1,ipiv,s0,dd+1,info)
682        do i=1,ngr
683         do j=1,nu
684          do j1=1,nvv
685           s1=0
686c          calculation sum for new solution
687           do jd=1,dd
688            s1=s1+s0(jd)*ggo(jd,j,j1,i)+lam*s0(jd)*dgg(jd,j,j1,i)
689           enddo
690c          change of  diis basis
691           do jd=1,dd-1
692            ggo(jd,j,j1,i)=ggo(jd+1,j,j1,i)
693            dgg(jd,j,j1,i)=dgg(jd+1,j,j1,i)
694           enddo
695           hnew(j,j1,i)=s1
696          enddo
697         enddo
698        enddo
699       endif
700c      matrix singular or has illegal value see dgeesv
701       if((k1.eq.dd).and.(info.ne.0)) then
702        k1=1
703        print*, info
704        go to 1
705       endif
706c      change count for diis
707       if(k1.lt.dd) then
708        k1=k1+1
709        else
710        k1=dd
711       endif
712 1     continue
713       return
714       end subroutine
715c
716c
717       subroutine thermo(nu,nv,nvv,ngr,icl,ims,den,t,rgrid,qqu,qvv,
718     * hnew,cnew,solperm,tau,muh,mugf)
719       implicit none
720#include "stdio.fh"
721       real*8 muh, mugf, muuh(1:nu,1:nvv),muugf(1:nu,1:nvv)
722       integer  nu, nv,nvv, ngr, i,j1,j,icl,nd
723       real*8 rgrid(1:ngr),dentv(1:nvv),cu(1:nu*nvv,1:ngr)
724       real*8 pi,dr,kb,t,kcal,theta,ck,echar,solperm,nav,tau
725       real*8 den(1:nv),gu(1:nu*nvv,1:ngr),nc(1:nu*nvv,1:ngr)
726       integer ims(1:nv)
727       real*8 cnew(1:nu,1:nvv,1:ngr), hnew(1:nu,1:nvv,1:ngr)
728       real*8 ul(1:nu,1:nvv,1:ngr),hq(1:nu,1:nvv,1:ngr)
729       real*8 ulr(1:nu,1:nvv,1:ngr),qvv(1:nvv),qqu(1:nu)
730       external util_erf
731       double precision util_erf
732       do j=1,nv
733        dentv(ims(j))=den(ims(j))
734       enddo
735       pi=2*asin(1.0)
736       kb=8.31441e-3
737       kcal=0.0019859
738       nd=(nvv+1)*nvv/2
739       dr=rgrid(2)-rgrid(1)
740       do j=1,nu
741        do j1=1,nvv
742         gu(nvv*(j-1)+j1,1)=0
743         nc(nvv*(j-1)+j1,1)=0
744         cu(nvv*(j-1)+j1,1)=cnew(j,j1,1)
745         do i=2,ngr
746          cu(nvv*(j-1)+j1,i)=cnew(j,j1,i)
747          gu(nvv*(j-1)+j1,i)=hnew(j,j1,i-1)+1
748          nc(nvv*(j-1)+j1,i)=nc(nvv*(j-1)+j1,i-1)
749     *    +4*pi*dr*dentv(j1)*gu(nvv*(j-1)+j1,i)*rgrid(i)**2.0
750         enddo
751        enddo
752       enddo
753       open(3, file='rdf_out.data', status='unknown')
754        do i=1,ngr
755         write(3,*) rgrid(i),(gu(j,i), j=1,nu*nvv)
756        enddo
757       close(3)
758       open(3, file='nc_out.data', status='unknown')
759        do i=1,ngr
760         write(3,*) rgrid(i),(nc(j,i), j=1,nu*nvv)
761        enddo
762       close(3)
763       open(3, file='c_out.data', status='unknown')
764        do i=1,ngr
765         write(3,*) rgrid(i),(cu(j,i), j=1,nu*nvv)
766        enddo
767       close(3)
768       mugf=0
769       muh=0
770       do j=1,nu
771        do j1=1,nvv
772         muugf(j,j1)=0
773         muuh(j,j1)=0
774        enddo
775       enddo
776       dr=rgrid(2)-rgrid(1)
777       nav=6.02214179e+23
778       echar=4.803e-10
779       ck=nav*echar**2/solperm*0.01
780c      ulr electrostatic part of the interaction potential in real space
781       do j=1,nu
782        do j1=1,nvv
783         do i=2,ngr
784          ulr(j,j1,i)=ck*qqu(j)*qvv(j1)/rgrid(i)*util_erf(tau*rgrid(i))
785         enddo
786         ulr(j,j1,1)=ulr(j,j1,2)
787        enddo
788       enddo
789c      calculations of partial contibutions
790       do j=1,nu
791        do j1=1,nvv
792         do i=2,ngr
793          if(icl.eq.1) then
794            hq(j,j1,i)=hnew(j,j1,i)
795          endif
796          if(icl.eq.2) then
797           hq(j,j1,i)= -theta(-hnew(j,j1,i))
798          endif
799          muugf(j,j1)=muugf(j,j1)-(2*cnew(j,j1,i)+hnew(j,j1,i)*
800     *    (cnew(j,j1,i)-ulr(j,j1,i)/kb/t))*rgrid(i)**2*dr*2*pi
801          muuh(j,j1)=muuh(j,j1)-(2*cnew(j,j1,i)+hnew(j,j1,i)*
802     *    (cnew(j,j1,i)-ulr(j,j1,i)/kb/t))*rgrid(i)**2*dr*2*pi
803     *    +hnew(j,j1,i)*hq(j,j1,i)*rgrid(i)**2*dr*2*pi
804         enddo
805        enddo
806       enddo
807       do j=1,nu
808        do j1=1,nv
809         mugf=mugf+den(j1)*muugf(j,ims(j1))*kcal*t
810         muh=muh+den(j1)*muuh(j,ims(j1))*kcal*t
811        enddo
812       enddo
813
814c      mugf chem.potential in gaussian approximation
815c      muh   chem.potential in hnc approximation
816       write(luout,108) "(hnc approximation)",muh
817       write(luout,108) "(gaussian approximation)",mugf
818  108  format("Chemical potential",A,T43,F10.4)
819       return
820       end subroutine
821
822c
823      real*8 function theta(x)
824      real*8 x
825       if(x.le.0.0)then
826        theta=x
827        else
828        theta=0
829       endif
830       return
831       end
832c
833
834c
835c uses realft
836c calculates the sine transform of a set of n real-valued data points stored in array y(1:n).
837c the number n must be a power of 2. on exit y is replaced by its transform. this program,
838c without changes, also calculates the inverse sine transform, but in this case the output array
839c should be multiplied by 2/n.
840c
841      subroutine sinft(y,n)
842      integer n
843      real*8 y(n)
844      integer j
845      real*8 sum,y1,y2
846      double precision theta,wi,wpi,wpr,wr,wtemp
847      theta=3.141592653589793d0/dble(n)
848      wr=1.0d0
849      wi=0.0d0
850      wpr=-2.0d0*sin(0.5d0*theta)**2
851      wpi=sin(theta)
852      y(1)=0.0
853      do 11 j=1,n/2
854        wtemp=wr
855        wr=wr*wpr-wi*wpi+wr
856        wi=wi*wpr+wtemp*wpi+wi
857        y1=wi*(y(j+1)+y(n-j+1))
858        y2=0.5*(y(j+1)-y(n-j+1))
859        y(j+1)=y1+y2
860        y(n-j+1)=y1-y2
86111    continue
862      call realft(y,n,+1)
863      sum=0.0
864      y(1)=0.5*y(1)
865      y(2)=0.0
866      do 12 j=1,n-1,2
867        sum=sum+y(j)
868        y(j)=y(j+1)
869        y(j+1)=sum
87012    continue
871      return
872      end subroutine
873c
874cc uses four1
875c calculates the fourier transform of a set of n real-valued data points. replaces this data
876c (which is stored in array data(1:n)) by the positive frequency half of its complex fourier
877c transform. the real-valued rst and last components of the complex transform are returned
878c as elements data(1) and data(2), respectively. n must be a power of 2. this routine
879c also calculates the inverse transform of a complex data array if it is the transform of real
880c data. (result in this case must be multiplied by 2/n.)
881c
882      subroutine realft(data,n,isign)
883      integer isign,n
884      real*8 data(n)
885      integer i,i1,i2,i3,i4,n2p3
886      real*8 c1,c2,h1i,h1r,h2i,h2r,wis,wrs
887      double precision theta,wi,wpi,wpr,wr,wtemp
888      theta=3.141592653589793d0/dble(n/2)
889      c1=0.5
890      if (isign.eq.1) then
891        c2=-0.5
892        call four1(data,n/2,+1)
893      else
894        c2=0.5
895        theta=-theta
896      endif
897      wpr=-2.0d0*sin(0.5d0*theta)**2
898      wpi=sin(theta)
899      wr=1.0d0+wpr
900      wi=wpi
901      n2p3=n+3
902      do 11 i=2,n/4
903        i1=2*i-1
904        i2=i1+1
905        i3=n2p3-i2
906        i4=i3+1
907        wrs=sngl(wr)
908        wis=sngl(wi)
909        h1r=c1*(data(i1)+data(i3))
910        h1i=c1*(data(i2)-data(i4))
911        h2r=-c2*(data(i2)+data(i4))
912        h2i=c2*(data(i1)-data(i3))
913        data(i1)=h1r+wrs*h2r-wis*h2i
914        data(i2)=h1i+wrs*h2i+wis*h2r
915        data(i3)=h1r-wrs*h2r+wis*h2i
916        data(i4)=-h1i+wrs*h2i+wis*h2r
917        wtemp=wr
918        wr=wr*wpr-wi*wpi+wr
919        wi=wi*wpr+wtemp*wpi+wi
92011    continue
921      if (isign.eq.1) then
922        h1r=data(1)
923        data(1)=h1r+data(2)
924        data(2)=h1r-data(2)
925      else
926        h1r=data(1)
927        data(1)=c1*(h1r+data(2))
928        data(2)=c1*(h1r-data(2))
929        call four1(data,n/2,-1)
930      endif
931      return
932      end subroutine
933cc replaces data(1:2*nn) by its discrete fourier transform, if isign is input as 1; or replaces
934c data(1:2*nn) by nn times its inverse discrete fourier transform, if isign is input as -1.
935c data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn
936c must be an integer power of 2 (this is not checked for!).
937      subroutine four1(data,nn,isign)
938      integer isign,nn
939      real*8 data(2*nn)
940      integer i,istep,j,m,mmax,n
941      real*8 tempi,tempr
942      double precision theta,wi,wpi,wpr,wr,wtemp
943      n=2*nn
944      j=1
945      do 11 i=1,n,2
946        if(j.gt.i)then
947          tempr=data(j)
948          tempi=data(j+1)
949          data(j)=data(i)
950          data(j+1)=data(i+1)
951          data(i)=tempr
952          data(i+1)=tempi
953        endif
954        m=n/2
9551       if ((m.ge.2).and.(j.gt.m)) then
956          j=j-m
957          m=m/2
958        goto 1
959        endif
960        j=j+m
96111    continue
962      mmax=2
9632     if (n.gt.mmax) then
964        istep=2*mmax
965        theta=6.28318530717959d0/(isign*mmax)
966        wpr=-2.d0*sin(0.5d0*theta)**2
967        wpi=sin(theta)
968        wr=1.d0
969        wi=0.d0
970        do 13 m=1,mmax,2
971          do 12 i=m,n,istep
972            j=i+mmax
973            tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
974            tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
975            data(j)=data(i)-tempr
976            data(j+1)=data(i+1)-tempi
977            data(i)=data(i)+tempr
978            data(i+1)=data(i+1)+tempi
97912        continue
980          wtemp=wr
981          wr=wr*wpr-wi*wpi+wr
982          wi=wi*wpr+wtemp*wpi+wi
98313      continue
984        mmax=istep
985      goto 2
986      endif
987      return
988      end subroutine
989
990c $Id$
991