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