1 subroutine prp_step(mdstep,stime,eww,esw,ess,fss,esk, 2 + epme,uqmd,uqmmm) 3c 4c $Id$ 5c 6 implicit none 7c 8#include "prp_common.fh" 9#include "msgids.fh" 10#include "mafdecls.fh" 11#include "global.fh" 12c 13 logical frequency 14 external frequency 15c 16 integer mdstep 17 real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),epme(3) 18 real*8 fss(msf,msf,3,2) 19 real*8 stime,esk(msf),uqmd,uqmmm 20c 21 if(me.eq.0) call prp_stat(mdstep,stime,eww,esw,ess,esk,epme(iset)) 22c 23 if(frequency(mdstep,npener)) then 24 call cf_print_energy(lfnout) 25 endif 26c 27 if(frequency(mdstep,nfprop)) call prp_record() 28c 29 return 30 end 31 subroutine prp_proper(mdstep,stime,eww,esw,ess,fss,esk,epme, 32 + uqmd,uqmmm,epot,epotw,epotsw,epots,vol,dwr,dsr,ekin,etot, 33 + npolit,gsm,esa,box,xsm) 34c 35 implicit none 36c 37#include "prp_common.fh" 38#include "msgids.fh" 39#include "mafdecls.fh" 40#include "global.fh" 41c 42 external timer_wall,timer_wall_total 43 real*8 timer_wall,timer_wall_total 44c 45 integer mdstep 46 real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),epme(3) 47 real*8 fss(msf,msf,3,2) 48 real*8 epot,epotw,epotsw,epots,vol,ekin,etot,epmec,ubias 49 real*8 eq,el,ep,ek,et,ewp,uqmd,uqmmm,box(3) 50 real*8 tempw,temps,stime,esk(msf),ewk,ep2,ep3,edrs,epmf 51 integer nwwl,nwws,nswl,nsws,nssl,nsss,nshitw,nshits,npolit 52 real*8 dwr,dsr(msm),gsm(msm,4,2),esa(nsa,2) 53 real*8 prest(3,3),virt(3,3),ekct(3,3),ep2m,ep3m,epmec2,epmec3 54 real*8 xsm(msm,3),dx(3),upmf(100) 55c 56 integer i,j,k,it,lenp 57c 58 call timer_start(54) 59c 60 if(ltwin) then 61 do 1 i=1,mpe 62 eww(i,1)=eww(i,1)+eww(i,2) 63 1 continue 64 do 2 i=1,msf 65 do 3 j=1,8 66 esw(i,j,1)=esw(i,j,1)+esw(i,j,2) 67 3 continue 68 2 continue 69 do 4 i=1,msf 70 do 5 j=1,msf 71 do 6 k=1,8 72 ess(i,j,k,1)=ess(i,j,k,1)+ess(i,j,k,2) 73 6 continue 74 fss(i,j,1,1)=fss(i,j,1,1)+fss(i,j,1,2) 75 fss(i,j,2,1)=fss(i,j,2,1)+fss(i,j,2,2) 76 fss(i,j,3,1)=fss(i,j,3,1)+fss(i,j,3,2) 77 5 continue 78 4 continue 79 endif 80c 81 el=eww(5,1)+eww(7,1) 82 eq=eww(6,1)+eww(8,1) 83 ep=zero 84 ewp=zero 85 do 9 i=1,mpe 86 ep=ep+eww(i,1) 87 ewp=ewp+eww(i,1) 88 do 10 j=1,msf 89 ep=ep+esw(j,i,1) 90 do 11 k=1,msf 91 ep=ep+ess(k,j,i,1) 92 11 continue 93 10 continue 94 9 continue 95 epmec=eww(9,1) 96 epmec2=eww(9,1) 97 epmec3=eww(9,1) 98 do 12 j=1,msf 99 el=el+esw(j,5,1) 100 eq=eq+esw(j,6,1) 101 do 13 k=1,msf 102 el=el+ess(k,j,5,1)+ess(k,j,7,1) 103 eq=eq+ess(k,j,6,1)+ess(k,j,8,1) 104 epmec=epmec+ess(k,j,9,1) 105 epmec2=epmec2+ess(k,j,10,1) 106 epmec3=epmec3+ess(k,j,11,1) 107 13 continue 108 12 continue 109c 110 if(me.eq.0) ep=ep+uqmd+uqmmm 111c 112c if using cafe get properties from it 113c 114 if(lcafe) then 115 call cf_proper(volume,temp,tempw,temps,pres,tmpscl,tmsscl,prsscl, 116 + ewk,nwwl,nwws,nswl,nsws,nssl,nsss,nshitw,nshits, 117 + ep2,ep3,ep2m,ep3m,edrs,epmf,virial,prest,virt,ekct,ubias,upmf) 118 endif 119c 120c if using space get properties from it 121c 122 if(lspac) then 123 endif 124c 125 if(me.eq.0) ep=ep+edrs 126 ep=ep+epmf 127c 128 if(lpme) then 129 ep=ep+epme(iset) 130 eq=eq+epme(iset) 131 ep2=ep2+epme(2)-epme(1)+epmec2-epmec 132 ep3=ep3+epme(3)-epme(1)+epmec3-epmec 133 endif 134c 135 ek=ewk 136 do 111 i=1,msf 137 ek=ek+esk(i) 138 111 continue 139 et=ep+ek 140c 141c fill the property vector 142c 143 do 7 i=1,maxpro 144 p(i)=zero 145 7 continue 146c 147 if(me.eq.0) then 148c 149 if(ntype.eq.0) call sp_gettp(temp,pres) 150c 151 p(1)=dble(mdstep) 152 p(2)=stime 153 endif 154 p(3)=dble(nwws) 155 p(4)=dble(nwwl) 156 p(5)=dble(nsws) 157 p(6)=dble(nswl) 158 p(7)=dble(nsss) 159 p(8)=dble(nssl) 160 p(9)=dble(nshitw) 161 p(10)=dble(nshits) 162 p(11)=dble(npolit) 163 p(33)=zero 164 p(34)=zero 165 if(me.eq.0) then 166 p(12)=volume 167 p(13)=1.6605655d0*wbox/volume 168 p(14)=pres 169 p(15)=prsscl 170 p(16)=temp 171 p(17)=tempw 172 p(18)=temps 173 p(21)=tmpscl 174 p(22)=tmsscl 175 p(27)=ek 176 p(32)=ek 177 p(33)=virial 178 p(34)=pres*volume 179 p(64)=ewk 180 p(67)=ewk 181 endif 182 p(24)=eq 183 p(25)=el 184 p(26)=ep 185 p(32)=p(32)+ep 186 p(36)=ep2 187 p(37)=ep3 188 p(38)=ep2+ep2m 189 p(39)=ep3+ep3m 190 if(me.eq.0) p(40)=dfree 191 p(50)=ubias 192 p(52)=eww(6,1) 193 p(53)=eww(5,1) 194 p(54)=eww(8,1) 195 p(55)=eww(7,1) 196 p(56)=eww(1,1) 197 p(57)=eww(2,1) 198 p(58)=eww(13,1) 199 p(59)=eww(3,1) 200 p(60)=eww(4,1) 201 p(67)=ewp 202 p(68)=p(68)+ewp 203 p(70)=dwr 204 if(lnoe) p(75)=edrs 205 if(lpmf) p(76)=epmf 206 if(ntype.ne.3) then 207 p(77)=epme(iset) 208 p(80)=epmec 209 else 210 p(77)=epme(1) 211 p(78)=epme(2)-p(77) 212 p(79)=epme(3)-p(77) 213 p(80)=epme(1)+epmec 214 p(81)=epme(2)+epmec2-p(80) 215 p(82)=epme(3)+epmec3-p(80) 216 endif 217 p(84)=timer_wall(202) 218 p(85)=timer_wall(203) 219 if(me.eq.0) p(86)=timer_wall(203) 220 if(me.eq.0) then 221 p(87)=virt(1,1) 222 p(88)=virt(1,2) 223 p(89)=virt(1,3) 224 p(90)=virt(2,1) 225 p(91)=virt(2,2) 226 p(92)=virt(2,3) 227 p(93)=virt(3,1) 228 p(94)=virt(3,2) 229 p(95)=virt(3,3) 230 p(96)=prest(1,1) 231 p(97)=prest(1,2) 232 p(98)=prest(1,3) 233 p(99)=prest(2,1) 234 p(100)=prest(2,2) 235 p(101)=prest(2,3) 236 p(102)=prest(3,1) 237 p(103)=prest(3,2) 238 p(104)=prest(3,3) 239 p(105)=ekct(1,1) 240 p(106)=ekct(1,2) 241 p(107)=ekct(1,3) 242 p(108)=ekct(2,1) 243 p(109)=ekct(2,2) 244 p(110)=ekct(2,3) 245 p(111)=ekct(3,1) 246 p(112)=ekct(3,2) 247 p(113)=ekct(3,3) 248 p(114)=box(1) 249 p(115)=box(2) 250 p(116)=box(3) 251 endif 252c 253 maxp=isprop 254 it=isprop 255c 256 if(nsf.gt.0) then 257 do 33 i=1,nsf 258 it=isprop+(i-1)*30 259 p(it+2)=ess(i,i,6,1) 260 p(it+3)=ess(i,i,5,1) 261 p(it+5)=esw(i,6,1) 262 p(it+6)=esw(i,5,1) 263 p(it+7)=ess(i,i,1,1) 264 p(it+8)=ess(i,i,2,1) 265 p(it+9)=ess(i,i,13,1) 266 p(it+10)=ess(i,i,3,1) 267 p(it+11)=ess(i,i,4,1) 268 p(it+14)=ess(i,i,7,1) 269 p(it+15)=ess(i,i,8,1) 270 p(it+16)=ess(i,i,5,1) 271 p(it+17)=ess(i,i,6,1) 272 if(me.eq.0) then 273 p(it+19)=esk(i) 274 endif 275 p(it+20)=ess(i,i,5,1)+ess(i,i,6,1)+ess(i,i,7,1)+ess(i,i,8,1)+ 276 + ess(i,i,1,1)+ess(i,i,2,1)+ 277 + ess(i,i,3,1)+ess(i,i,4,1)+half*(esw(i,6,1)+esw(i,5,1)) 278 do 34 j=i+1,nsf 279 p(it+20)=p(it+18)+half*(ess(i,j,6,1)+ess(j,i,6,1)+ 280 + ess(i,j,5,1)+ess(j,i,5,1)) 281 34 continue 282 33 continue 283 it=it+30 284 if(nsf.gt.1) then 285 do 35 i=1,nsf-1 286 do 36 j=i+1,nsf 287 it=it+1 288 it=it+1 289 p(it)=ess(i,j,6,1)+ess(j,i,6,1) 290 it=it+1 291 p(it)=ess(i,j,5,1)+ess(j,i,5,1) 292 it=it+1 293 it=it+1 294 p(it)=fss(i,j,1,1)-fss(j,i,1,1) 295 it=it+1 296 p(it)=fss(i,j,2,1)-fss(j,i,2,1) 297 it=it+1 298 p(it)=fss(i,j,3,1)-fss(j,i,3,1) 299 it=it+1 300 p(it)=sqrt((fss(i,j,1,1)-fss(j,i,1,1))**2+ 301 + (fss(i,j,2,1)-fss(j,i,2,1))**2+(fss(i,j,3,1)-fss(j,i,3,1))**2) 302 36 continue 303 35 continue 304 endif 305 maxp=it 306c 307 if(nsm.gt.0.and.npstat.gt.0) then 308 it=maxp 309 do 37 i=1,nsm 310 it=it+1 311 if(me.eq.0) p(it)=gsm(i,iset,1) 312 if(p(it).lt.tiny) lp(it)=.false. 313 37 continue 314 do 55 i=1,nsm 315 it=it+1 316 if(me.eq.0) p(it)=gsm(i,4,1) 317 if(p(it).lt.tiny) lp(it)=.false. 318 55 continue 319 do 155 i=1,nsm-1 320 do 156 j=i+1,nsm 321 dx(1)=xsm(i,1)-xsm(j,1) 322 dx(2)=xsm(i,2)-xsm(j,2) 323 dx(3)=xsm(i,3)-xsm(j,3) 324 if(lpbc) call cf_pbc(1,dx,1,dx,1,0,1,1) 325 it=it+1 326 if(me.eq.0) p(it)=sqrt(dx(1)**2+dx(2)**2+dx(3)**2) 327 156 continue 328 155 continue 329 endif 330c 331 endif 332c 333 if(npmfi.gt.0) then 334 do 356 i=1,npmfi 335 it=it+1 336 p(it)=upmf(i) 337 356 continue 338 endif 339c 340 call timer_stop(54) 341c 342 if(iprof.eq.1) then 343 do 56 i=1,55 344 it=it+1 345 p(it)=timer_wall_total(i) 346 56 continue 347 endif 348c 349 maxp=it 350c 351 if(npener.gt.0) call cf_add_esa(esa) 352 if(np.gt.0) then 353 call ga_dgop(mrg_d44,p(9),3,'max') 354 if(me.ne.0) then 355 p(9)=zero 356 p(10)=zero 357 p(11)=zero 358 endif 359 if(maxp+mpe*(1+msf*(1+msf)).gt.maxpro) then 360 call ga_dgop(mrg_d45,p,maxp,'+') 361 call ga_dgop(mrg_d40,eww,mpe,'+') 362 call ga_dgop(mrg_d41,esw,msf*mpe,'+') 363 call ga_dgop(mrg_d42,ess,msf*msf*mpe,'+') 364 else 365 lenp=maxp 366 do 501 i=1,mpe 367 p(lenp+i)=eww(i,1) 368 501 continue 369 lenp=lenp+mpe 370 do 502 i=1,mpe 371 do 503 j=1,msf 372 lenp=lenp+1 373 p(lenp)=esw(j,i,1) 374 503 continue 375 502 continue 376 do 504 i=1,mpe 377 do 505 j=1,msf 378 do 506 k=1,msf 379 lenp=lenp+1 380 p(lenp)=ess(k,j,i,1) 381 506 continue 382 505 continue 383 504 continue 384 call ga_dgop(mrg_d45,p,lenp,'+') 385 lenp=maxp 386 do 511 i=1,mpe 387 eww(i,1)=p(lenp+i) 388 511 continue 389 lenp=lenp+mpe 390 do 512 i=1,mpe 391 do 513 j=1,msf 392 lenp=lenp+1 393 esw(j,i,1)=p(lenp) 394 513 continue 395 512 continue 396 do 514 i=1,mpe 397 do 515 j=1,msf 398 do 516 k=1,msf 399 lenp=lenp+1 400 ess(k,j,i,1)=p(lenp) 401 516 continue 402 515 continue 403 514 continue 404 endif 405cxxxxx call ga_dgop(mrg_d43,epme,3,'+') 406 if(npener.gt.0) call ga_dgop(mrg_d47,esa,2*nsa,'+') 407 endif 408c 409 if(p(85).ne.0.0d0) then 410 p(85)=(p(85)-p(84))/p(85) 411 endif 412c 413 if(temp.eq.zero.or.me.ne.0) then 414 p(117)=zero 415 p(118)=zero 416 p(119)=zero 417 p(120)=zero 418 p(51)=zero 419 else 420 p(117)=exp(-p(36)/(rgas*temp)) 421 p(118)=exp(-p(37)/(rgas*temp)) 422 p(119)=exp(-p(38)/(rgas*temp)) 423 p(120)=exp(-p(39)/(rgas*temp)) 424 p(51)=exp(-p(50)/(rgas*temp)) 425 if(nbias.gt.0) p(41)=p(40)*p(51) 426 endif 427c 428 if(nwm.gt.0) then 429 dwr=p(70)/dble(nwm) 430 p(70)=dwr 431 if(stime.ne.0) then 432 p(71)=1.0d-6*dwr/(6.0d0*stime) 433 else 434 p(71)=zero 435 endif 436 endif 437c 438 epot=p(26) 439 epots=zero 440 epotsw=zero 441 epotw=eww(5,1)+eww(6,1)+eww(8,1) 442 do 40 i=1,msf 443 epotsw=epotsw+esw(i,5,1)+esw(i,6,1)+esw(i,8,1) 444 do 41 j=1,msf 445 epots=epots+ess(i,j,1,1)+ess(i,j,2,1)+ess(i,j,3,1)+ess(i,j,4,1) 446 epots=epots+ess(i,j,5,1)+ess(i,j,6,1)+ess(i,j,8,1) 447 41 continue 448 40 continue 449c 450 epots=zero 451 epotsw=zero 452 epotw=zero 453 do 42 i=1,mpe 454 epotw=epotw+eww(i,1) 455 do 38 j=1,msf 456 epotsw=epotsw+esw(j,i,1) 457 do 39 k=1,msf 458 epots=epots+ess(k,j,i,1) 459 39 continue 460 38 continue 461 42 continue 462c 463 vol=volume 464 ekin=p(27) 465 etot=ekin+epot 466c 467 return 468 end 469 subroutine prp_stat(mdstep,stime,eww,esw,ess,esk,epme) 470c 471 implicit none 472c 473#include "prp_common.fh" 474#include "msgids.fh" 475#include "mafdecls.fh" 476#include "global.fh" 477c 478 logical frequency 479 external frequency 480c 481 integer mdstep 482 real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),epme 483 real*8 stime,esk(msf) 484c 485 integer i,j 486 character*10 pdate,ptime 487 real*8 facs,fact,tfacs,tfact,rt 488c 489 nsum=nsum+1 490 nsumt=nsumt+1 491 nsump=nsump+1 492 do 8 i=1,maxpro 493 if(abs(p(i)).lt.tiny) p(i)=zero 494 psum(i)=psum(i)+p(i) 495 p2sum(i)=p2sum(i)+p(i)*p(i) 496 pslop(i)=pslop(i)+stime*p(i) 497 psumt(i)=psumt(i)+p(i) 498 p2sumt(i)=p2sumt(i)+p(i)*p(i) 499 pslopt(i)=pslopt(i)+stime*p(i) 500 psump(i)=psump(i)+p(i) 501 8 continue 502 tsum=tsum+stime 503 t2sum=t2sum+stime*stime 504 tsumt=tsumt+stime 505 t2sumt=t2sumt+stime*stime 506c 507 if(frequency(mdstep,nfoutp)) then 508 if(.not.lhdr) then 509 call swatch(pdate,ptime) 510 write(lfnout,1000) pdate,ptime 511 1000 format(/,' MOLECULAR DYNAMICS TIME STEP INFORMATION',T110,2A10,//, 512 + ' Time Temp Pres Volume Tscalw Tscals Pscal ', 513 + ' U(ele) U(vdW) U(pot) U(kin) U(tot) ',/, 514 + ' ps K Pa nm**3 ', 515 + ' kJ/mol kJ/mol kJ/mol kJ/mol kJ/mol ',/) 516 if(lpstep) write(lfnout,1010) 517 1010 format(14X, 518 + ' U(bnd) U(ang) U(dih) U(imp) ', 519 + ' Ui(3rd) Ui(non) ', 520 + ' Uw(ele) Uw(vdW) U(kin) U(pot) U(tot)',/,14X, 521 + ' kJ/mol kJ/mol kJ/mol kJ/mol ', 522 + ' kJ/mol kJ/mol ', 523 + ' kJ/mol kJ/mol kJ/mol kJ/mol kJ/mol ',/) 524 lhdr=.true. 525 endif 526 write(lfnout,1001) stime,temp,pres,volume,p(21),p(22),p(15), 527 + p(24),p(25),p(26),p(27),p(32) 528 1001 format(1x,f10.5,0pf8.2,1pe9.2,0pf10.3,3f7.4,5(1pe11.4),i5,i7) 529 if(lpstep) then 530 if(nwm.gt.0) then 531 rt=one/nwm 532 write(lfnout,1011) rt*p(56),rt*p(57),rt*p(58),rt*p(59),rt*p(52), 533 + rt*p(53),rt*p(54),rt*p(55),rt*p(64),rt*p(66),rt*p(67) 534 1011 format(' solvent ',11f11.2) 535 endif 536 do 1 i=1,nsf 537 j=isprop+(i-1)*30 538 rt=p(j+7)+p(j+8)+p(j+9)+p(j+10)+p(j+2)+p(j+3)+p(j+5)+p(j+6) 539 write(lfnout,1012) i,p(j+7),p(j+8),p(j+9),p(j+10), 540 + p(j+14)+p(j+15),p(j+2)+p(j+3),p(j+5),p(j+6),p(j+19),rt,rt+p(j+19) 541 1012 format(' solute',i3,11f11.2) 542 1 continue 543 endif 544 endif 545c 546 if(frequency(mdstep,nfstat)) then 547 call swatch(pdate,ptime) 548 write(lfnout,2000) pdate,ptime,nsum,nsumt 549 2000 format(/,' MOLECULAR DYNAMICS STATISTICAL INFORMATION',t110,2a10, 550 + //,t41,2(3X,'Statistics over last ',I8,' steps',2X),/, 551 + t41,2(3X,'Average',5X,'RMS fluct',5X,'Drift/ps',3X),/) 552 facs=one/dble(nsum) 553 fact=one/dble(nsumt) 554 tfacs=one/(t2sum-facs*tsum*tsum) 555 tfact=one/(t2sumt-fact*tsumt*tsumt) 556 do 15 i=1,nprop 557 j=ixp(i) 558 if(lp(j)) write(lfnout,2001) pronam(j)(1:39), 559 + psum(j)*facs,sqrt(abs((p2sum(j)-psum(j)*psum(j)*facs)*facs)), 560 + (pslop(j)-facs*psum(j)*tsum)*tfacs, 561 + psumt(j)*fact,sqrt(abs((p2sumt(j)-psumt(j)*psumt(j)*fact)*fact)), 562 + (pslopt(j)-fact*psumt(j)*tsumt)*tfact, 563 + pronam(j)(40:50) 564 2001 format(1x,a39,t41,2(3(1pe12.5,1x),1x),a11) 565 psum(j)=zero 566 p2sum(j)=zero 567 pslop(j)=zero 568 15 continue 569 tsum=zero 570 t2sum=zero 571 nsum=0 572 lhdr=.false. 573 endif 574c 575 return 576 end 577 subroutine prp_print() 578c 579 implicit none 580c 581#include "prp_common.fh" 582c 583 character*10 pdate,ptime 584 integer i,j 585c 586 if(me.ne.0) return 587c 588 call swatch(pdate,ptime) 589 write(lfnout,1000) pdate,ptime 590 1000 format(/,' SINGLE POINT PROPERTIES',t110,2a10,/) 591 do 1 i=1,nprop 592 j=ixp(i) 593 if(lp(j)) write(lfnout,1001) pronam(j)(1:39),p(j),pronam(j)(40:50) 594 1001 format(1x,a39,t41,1pe18.9,1x,a11) 595 1 continue 596c 597 return 598 end 599 subroutine prp_record() 600c 601 implicit none 602c 603#include "prp_common.fh" 604c 605 integer i 606 character*10 pdate,ptime 607c 608 if(me.ne.0) return 609c 610 if(.not.lfhdr) then 611 call swatch(pdate,ptime) 612 write(lfnprp,3000) nprop,pdate,ptime,np,npfft 613 3000 format(i7,1x,2a10,2i5) 614 write(lfnprp,3001) (pronam(ixp(i)),i=1,nprop) 615 3001 format(a50) 616 lfhdr=.true. 617 endif 618 write(lfnprp,3002) 619 3002 format('frame') 620 if(iprop.eq.0) then 621 write(lfnprp,3003) (p(ixp(i)),i=1,nprop) 622 else 623 write(lfnprp,3003) (psump(ixp(i))/dble(nsump),i=1,nprop) 624 3003 format(4(1pe12.5)) 625 endif 626 nsump=0 627 do 1 i=1,nprop 628 psump(ixp(i))=zero 629 1 continue 630c 631 return 632 end 633 subroutine prp_header() 634c 635 implicit none 636c 637#include "prp_common.fh" 638c 639 lfhdr=.false. 640c 641 return 642 end 643 logical function prp_mcti_step(ida,lda) 644c 645 implicit none 646c 647#include "prp_common.fh" 648#include "mafdecls.fh" 649#include "msgids.fh" 650#include "global.fh" 651c 652 logical prp_mcti_acc 653 external prp_mcti_acc 654c 655 integer ida,lda 656c 657 real*8 fdata(28) 658 real*8 aver,drift,stderr,corerr,ratio 659 logical done 660c 661 lerror=ida.gt.lda 662c 663 call cf_mcti(fdata) 664c 665 call ga_dgop(mrg_d44,fdata,28,'+') 666c 667 if(me.eq.0) then 668 done=prp_mcti_acc(ida,dbl_mb(i_dfr),dbl_mb(i_dfrm),fdata, 669 + aver,drift,stderr,corerr,ratio) 670 endif 671c 672 if(np.gt.1) then 673 call ga_brdcst(mrg_d46,done,ma_sizeof(mt_log,1,mt_byte),0) 674 endif 675c 676 prp_mcti_step=done 677c 678 return 679 end 680 logical function prp_mcti_acc(ida,dfr,dfrm,fdata, 681 + aver,drift,stderr,corerr,ratio) 682c 683 implicit none 684c 685#include "prp_common.fh" 686c 687 integer ida 688 real*8 dfr(mda),dfrm(mda),fdata(28) 689 real*8 aver,drift,stderr,corerr,ratio,cerror 690c 691 integer i 692 logical done 693 real*8 dfrnom 694c 695 dfree=zero 696 dfrnom=zero 697 do 1 i=1,24 698 dfree=dfree+fdata(i) 699 if(i.ne.1.and.i.ne.13) dfrnom=dfrnom+fdata(i) 700 deriv(i)=deriv(i)+fdata(i) 701 1 continue 702 nderiv=nderiv+1 703c 704 dfr(ida)=dfree 705 dfrm(ida)=dfrnom 706 if(dfree.ne.dfrnom) lfreem=.true. 707 nda=ida 708c 709 2 continue 710c 711 if(lerror) then 712 call error(lauto,lappr,1000,dfr,ida, 713 + aver,drift,stderr,corerr,ratio) 714 cerror=corerr 715 if(.not.lauto) cerror=samrat*corerr 716 if(.not.lauto.and.cerror.lt.edacq) then 717 lauto=.true. 718 lappr=.true. 719 goto 2 720 else 721 if(lauto) samrat=ratio 722 lauto=.false. 723 lappr=.false. 724 endif 725 done=cerror.lt.edacq.and.drift.lt.ddacq 726 else 727 done=.false. 728 endif 729c 730 prp_mcti_acc=done 731 return 732 end 733 subroutine prp_mcti_run(rlambd,dlambd,ndec) 734c 735 implicit none 736c 737#include "prp_common.fh" 738#include "mafdecls.fh" 739c 740 real*8 rlambd,dlambd 741 integer ndec 742c 743 call prp_mcti_r(rlambd,dlambd,ndec,dbl_mb(i_dfr),dbl_mb(i_dfrm), 744 + psumt(16),psumt(117),psumt(118),psumt(119),psumt(120), 745 + psumt(51),psumt(41)) 746c 747 return 748 end 749 subroutine prp_mcti_r(rlambd,dlambd,ndec,dfr,dfrm,taver, 750 + ep2ave,ep3ave,ep2avm,ep3avm,ebias,dfbias) 751c 752 implicit none 753c 754#include "prp_common.fh" 755c 756 real*8 rlambd,dlambd,dfr(mda),dfrm(mda) 757 real*8 taver,ep2ave,ep3ave,ebias,dfbias,ep2avm,ep3avm 758 integer i,ndec 759c 760 if(me.eq.0) then 761 write(lfngib,1000) nderiv,nda,rlambd,dlambd,ndec,nsa, 762 + nbias,ebias/dble(nsumt),lfreem 763 1000 format(2i7,2f12.6,2i8,i4,e20.12,4x,l1) 764 write(lfngib,1001) deriv 765 1001 format(4e20.12) 766 write(lfngib,1002) (dfr(i),i=1,nda) 767 write(lfngib,1002) (dfrm(i),i=1,nda) 768 1002 format(4e20.12) 769 write(lfngib,1003) nsumt,taver/dble(nsumt), 770 + ep2ave/dble(nsumt),ep3ave/dble(nsumt),dfbias/ebias, 771 + ep2avm/dble(nsumt),ep3avm/dble(nsumt) 772 1003 format(i10,/,4e20.12,/,2e20.12) 773 endif 774 if(ndec.gt.0) call cf_wrtgib(lfngib) 775c 776 return 777 end 778 subroutine prp_mcti(n,filnam) 779c 780 implicit none 781c 782#include "prp_common.fh" 783#include "mafdecls.fh" 784#include "global.fh" 785c 786 integer n,i_dec,l_dec,ibl,i,igp_handle 787 character*80 filnam,string 788 character*5 gid 789c 790 if(npg.gt.1) then 791 if(me.eq.0) then 792 close(unit=lfngib,status='keep') 793 endif 794c 795c now go global 796c 797 igp_handle=ga_pgroup_get_default() 798 call ga_pgroup_set_default(ga_pgroup_get_world()) 799 call ga_sync() 800 if(ga_nodeid().eq.0) then 801c 802 ibl=index(filnam,' ')-1 803 string=filnam(1:ibl)//'.gib' 804 open(unit=lfngib,file=string,status='unknown',form='formatted') 805 rewind(lfngib) 806 do 1 i=1,npg 807 write(gid,3300) i-1 808 3300 format(i5.5) 809 string=filnam(1:ibl)//gid//'.gib' 810 open(unit=75,file=string,status='old',form='formatted') 811 rewind(unit=75) 812 2 continue 813 read(75,3301,end=99,err=99) string 814 write(lfngib,3301) string 815 3301 format(a) 816 goto 2 817 99 continue 818 close(unit=75) 819 1 continue 820 close(unit=lfngib,status='keep') 821c 822 endif 823c 824c go back to groups 825c 826 call ga_sync() 827 call ga_pgroup_set_default(igp_handle) 828c 829 else if(ipg.eq.npg-1.or.npg.eq.1) then 830 if(me.eq.0) then 831 close(unit=lfngib,status='keep') 832 endif 833 endif 834c 835 if(ipg.eq.npg-1.or.npg.eq.1) then 836c 837 if(me.eq.0) then 838 ibl=index(filnam,' ')-1 839 string=filnam(1:ibl)//'.gib' 840 open(unit=lfngib,file=string,status='old',form='formatted') 841 rewind(lfngib) 842 endif 843c 844 npgdec=n 845c 846 if(npgdec.gt.0) then 847 if(.not.ma_push_get(mt_dbl,6*nsa,'dec',l_dec,i_dec)) 848 + call md_abort('Failed to allocate dec',0) 849 else 850 if(.not.ma_push_get(mt_dbl,1,'dec',l_dec,i_dec)) 851 + call md_abort('Failed to allocate dec',0) 852 endif 853c 854 call prp_mcti_s(dbl_mb(i_dfr),dbl_mb(i_dfrm),dbl_mb(i_dec)) 855c 856 if(.not.ma_pop_stack(l_dec)) 857 + call md_abort('Failed to deallocate dec',0) 858c 859 endif 860c endif 861c 862 return 863 end 864 subroutine prp_mcti_s(dfr,dfrm,dec) 865c 866 implicit none 867c 868#include "prp_common.fh" 869#include "mafdecls.fh" 870c 871 real*8 dfr(mda),dfrm(mda),dec(6,nsa) 872c 873 character*10 pdate,ptime 874 integer i,j,k,number,ndec 875 real*8 rlambd,dlambd,ddrft,dsterr,dcerr,ratio 876 real*8 freeti,errti,drftti,taver,ep2ave,ep3ave,ebias,dfbias 877 real*8 ep2avm,ep3avm,dtmp 878 real*8 epr,epf,etp,slambd,rnum,fterm(24),gbias,freeb,freem 879 logical lfrm 880c 881 gbias=zero 882c 883 if(me.eq.0) then 884c 885 do 1 k=1,8 886c 887 if(k.le.4.and.npgdec.eq.0) goto 1 888c 889 call swatch(pdate,ptime) 890 if(k.eq.1) then 891 write(lfnout,1001) 892 1001 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ', 893 + 'DECOMPOSITION',//,' Solvent derivatives in kJ/mol',//, 894 + ' Run Lambda', 895 + ' Mass Solvnt LJ Solute LJ Solvnt el Solute el Bonds', 896 + ' Constrts Angles Dihedrals Impropers Slvnt pol Solut pol',/) 897 elseif(k.eq.2) then 898 write(lfnout,1002) 899 1002 format(//,' Solute derivatives in kJ/mol',//, 900 + ' Run Lambda', 901 + ' Mass Solvnt LJ Solute LJ Solvnt el Solute el Bonds', 902 + ' Constrts Angles Dihedrals Impropers Slvnt pol Solut pol',/) 903 elseif(k.eq.3) then 904 write(lfnout,1003) 905 1003 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ', 906 + 'DECOMPOSITION',//,' Solvent contributions in kJ/mol',//, 907 + ' Run Lambda', 908 + ' Mass Solvnt LJ Solute LJ Solvnt el Solute el Bonds', 909 + ' Constrts Angles Dihedrals Impropers Slvnt pol Solut pol',/) 910 elseif(k.eq.4) then 911 write(lfnout,1004) 912 1004 format(//,' Solute contributions in kJ/mol',//, 913 + ' Run Lambda', 914 + ' Mass Solvnt LJ Solute LJ Solvnt el Solute el Bonds', 915 + ' Constrts Angles Dihedrals Impropers Slvnt pol Solut pol',/) 916 elseif(k.eq.5) then 917 if(lfreem) then 918 write(lfnout,1005) pdate,ptime 919 1005 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ', 920 + 'EXCLUDING MASS CONTRIBUTIONS',t110,2a10,//, 921 + ' Run Lambda Size', 922 + ' Derivative Derivative Derivative Lambda', 923 + ' Free Energy Free Energy Free Energy Free Energy Sampling',/, 924 + 18x, 925 + ' Average Error Drift ', 926 + ' Accumulated Error Drift Corrected Ratio',/, 927 + 18x, 928 + ' kJ/mol kJ/mol kJ/mol ps ', 929 + ' kJ/mol kJ/mol kJ/mol ps kJ/mol',/) 930 else 931 write(lfnout,1006) pdate,ptime 932 1006 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION', 933 + t110,2a10,//, 934 + ' Run Lambda Size', 935 + ' Derivative Derivative Derivative Lambda', 936 + ' Free Energy Free Energy Free Energy Free Energy Sampling',/, 937 + 18x, 938 + ' Average Error Drift ', 939 + ' Accumulated Error Drift Corrected Ratio',/, 940 + 18x, 941 + ' kJ/mol kJ/mol kJ/mol ps ', 942 + ' kJ/mol kJ/mol kJ/mol ps kJ/mol',/) 943 endif 944 elseif(k.eq.6) then 945 if(lfreem) then 946 write(lfnout,1007) pdate,ptime 947 1007 format(//,' MULTICONFIGURATION THERMODYNAMIC INTEGRATION ', 948 + 'INCLUDING MASS CONTRIBUTIONS',t110,2a10,//, 949 + ' Run Lambda Size', 950 + ' Derivative Derivative Derivative Lambda', 951 + ' Free Energy Free Energy Free Energy Free Energy Sampling',/, 952 + 18x, 953 + ' Average Error Drift ', 954 + ' Accumulated Error Drift Corrected Ratio',/, 955 + 18x, 956 + ' kJ/mol kJ/mol kJ/mol ps ', 957 + ' kJ/mol kJ/mol kJ/mol ps kJ/mol',/) 958 endif 959 elseif(k.eq.7) then 960 if(lfreem) then 961 write(lfnout,1008) pdate,ptime 962 1008 format(//' MULTISTEP THERMODYNAMIC PERTURBATION ', 963 + 'EXCLUDING MASS CONTRIBUTIONS',t110,2a10,//, 964 + ' Run Lambda Size', 965 + ' Temperature Reverse Forward Ensemble Lambda', 966 + ' Accumulated Bias Corrected',/,18x, 967 + ' K kJ/mol kJ/mol kJ/mol ', 968 + ' kJ/mol kJ/mol kJ/mol',/) 969 else 970 write(lfnout,1009) pdate,ptime 971 1009 format(//' MULTISTEP THERMODYNAMIC PERTURBATION', 972 + t110,2a10,//, 973 + ' Run Lambda Size', 974 + ' Temperature Reverse Forward Ensemble Lambda', 975 + ' Accumulated Bias Corrected',/,18x, 976 + ' K kJ/mol kJ/mol kJ/mol ', 977 + ' kJ/mol kJ/mol kJ/mol',/) 978 endif 979 else 980 if(lfreem) then 981 write(lfnout,1010) pdate,ptime 982 1010 format(//' MULTISTEP THERMODYNAMIC PERTURBATION ', 983 + 'INCLUDING MASS CONTRIBUTIONS',t110,2a10,//, 984 + ' Run Lambda Size', 985 + ' Temperature Reverse Forward Ensemble Lambda', 986 + ' Accumulated Bias Corrected',/,18x, 987 + ' K kJ/mol kJ/mol kJ/mol ', 988 + ' kJ/mol kJ/mol kJ/mol',/) 989 endif 990 endif 991c 992 rewind(lfngib) 993c 994 freeti=zero 995 freeb=zero 996 freem=zero 997 errti=zero 998 drftti=zero 999 etp=zero 1000 slambd=zero 1001 gbias=zero 1002 do 2 i=1,24 1003 fterm(i)=zero 1004 2 continue 1005 if(npgdec.gt.0) then 1006 do 22 i=1,nsa 1007 dec(1,i)=zero 1008 dec(2,i)=zero 1009 dec(3,i)=zero 1010 dec(4,i)=zero 1011 dec(5,i)=zero 1012 dec(6,i)=zero 1013 22 continue 1014 endif 1015 do 3 i=1,mrun 1016 read(lfngib,2000) nderiv,nda,rlambd,dlambd,ndec,nsa,nbias,ebias, 1017 + lfrm 1018 if(lfrm) lfreem=.true. 1019 2000 format(2i7,2f12.6,2i8,i4,e20.12,4x,l1) 1020 read(lfngib,2001) deriv 1021 2001 format(4e20.12) 1022 if(mda.lt.nda) then 1023 read(lfngib,2002) (dtmp,j=1,nda) 1024 read(lfngib,2002) (dtmp,j=1,nda) 1025 else 1026 read(lfngib,2002) (dfr(j),j=1,nda) 1027 read(lfngib,2002) (dfrm(j),j=1,nda) 1028 endif 1029 2002 format(4e20.12) 1030 read(lfngib,2003) number,taver,ep2ave,ep3ave,dfbias, 1031 + ep2avm,ep3avm 1032 2003 format(i10,/,4e20.12,/,2e20.12) 1033 if(nbias.gt.0) gbias=-rgas*taver*log(ebias) 1034 if(ndec.gt.0) call cf_rdgib(lfngib,dec,dlambd/dble(ndec)) 1035c 1036 slambd=slambd+dlambd 1037 rnum=dlambd/dble(nderiv) 1038c 1039 do 4 j=1,24 1040 fterm(j)=fterm(j)+rnum*deriv(j) 1041 4 continue 1042c 1043 freem=freem+rnum*(deriv(1)+deriv(13)) 1044c 1045 if(k.eq.1) then 1046 rnum=one/dble(nderiv) 1047 write(lfnout,1101) i,slambd,(rnum*deriv(j),j=1,12) 1048 1101 format(i4,f7.3,12f10.3) 1049 elseif(k.eq.2) then 1050 rnum=one/dble(nderiv) 1051 write(lfnout,1102) i,slambd,(rnum*deriv(j),j=13,24) 1052 1102 format(i4,f7.3,12f10.3) 1053 elseif(k.eq.3) then 1054 write(lfnout,1103) i,slambd,(fterm(j),j=1,12) 1055 1103 format(i4,f7.3,12f10.5) 1056 elseif(k.eq.4) then 1057 write(lfnout,1104) i,slambd,(fterm(j),j=13,24) 1058 1104 format(i4,f7.3,12f10.5) 1059 elseif(k.eq.5) then 1060 call error(.true.,.true.,1000,dfrm,nda, 1061 + dfree,ddrft,dsterr,dcerr,ratio) 1062c 1063 freeti=freeti+dfree*dlambd 1064 freeb=freeb+dfbias*dlambd 1065 errti=errti+(dcerr*dlambd)**2 1066 drftti=drftti+ddrft*dlambd 1067c 1068 write(lfnout,1105) i,rlambd,nderiv,dfree,dcerr,ddrft/tstep, 1069 + slambd,freeti,sqrt(errti),drftti/tstep,freeb,ratio 1070 1105 format(i4,f7.3,i7,f12.3,2f12.6,f7.3,f12.3,2f12.6,f12.3,f10.5) 1071c 1072 elseif(k.eq.6) then 1073 if(lfreem) then 1074 call error(.true.,.true.,1000,dfr,nda, 1075 + dfree,ddrft,dsterr,dcerr,ratio) 1076c 1077 freeti=freeti+dfree*dlambd 1078 freeb=freeb+dfbias*dlambd 1079 errti=errti+(dcerr*dlambd)**2 1080 drftti=drftti+ddrft*dlambd 1081c 1082 write(lfnout,1106) i,rlambd,nderiv,dfree,dcerr,ddrft/tstep, 1083 + slambd,freeti,sqrt(errti),drftti/tstep,freeb,ratio 1084 1106 format(i4,f7.3,i7,f12.3,2f12.6,f7.3,f12.3,2f12.6,f12.3,f10.5) 1085 endif 1086c 1087 elseif(k.eq.7) then 1088 epr=-rgas*taver*log(ep2ave) 1089 epf=-rgas*taver*log(ep3ave) 1090 etp=etp-epr+epf 1091 write(lfnout,1107) i,rlambd,number,taver,epr,epf,epf-epr, 1092 + slambd,etp,-gbias,etp-gbias 1093 1107 format(i4,f7.3,i7,4f12.3,f7.3,3f12.3) 1094 else 1095 if(lfreem) then 1096 epr=-rgas*taver*log(ep2avm) 1097 epf=-rgas*taver*log(ep3avm) 1098 etp=etp-epr+epf 1099 write(lfnout,1108) i,rlambd,number,taver,epr,epf,epf-epr, 1100 + slambd,etp,-gbias,etp-gbias 1101 1108 format(i4,f7.3,i7,4f12.3,f7.3,3f12.3) 1102 endif 1103 endif 1104c 1105 3 continue 1106c 1107 if(ndec.gt.0.and.k.eq.4) call cf_print_deco(lfnout,dec) 1108c 1109 1 continue 1110c 1111 close(unit=lfngib) 1112c 1113 endif 1114c 1115 return 1116 end 1117 subroutine prp_wrtmro(lfnmro,ndec) 1118c 1119 implicit none 1120c 1121#include "prp_common.fh" 1122#include "mafdecls.fh" 1123c 1124 integer lfnmro,ndec 1125c 1126 call prp_wtmro(lfnmro,ndec,dbl_mb(i_dfr),dbl_mb(i_dfrm)) 1127c 1128 return 1129 end 1130 logical function prp_rdmri(lfnmri,ndec,mropt) 1131c 1132 implicit none 1133c 1134#include "prp_common.fh" 1135#include "mafdecls.fh" 1136c 1137 logical prp_rmri 1138 external prp_rmri 1139c 1140 integer lfnmri,ndec,mropt 1141c 1142 prp_rdmri=prp_rmri(lfnmri,ndec,mropt,dbl_mb(i_dfr),dbl_mb(i_dfrm)) 1143c 1144 return 1145 end 1146 subroutine prp_wtmro(lfnmro,ndec,dfr,dfrm) 1147c 1148 implicit none 1149c 1150#include "prp_common.fh" 1151c 1152 integer lfnmro,ndec 1153 real*8 dfr(mda),dfrm(mda) 1154c 1155 integer i 1156c 1157 if(me.eq.0) then 1158 write(lfnmro) nderiv,nda,nprop,nsum,nsumt,ndec,maxp 1159 write(lfnmro) deriv 1160 write(lfnmro) (dfr(i),i=1,nda) 1161 write(lfnmro) (dfrm(i),i=1,nda) 1162 write(lfnmro) tsum,t2sum,tsumt,t2sumt 1163 write(lfnmro) (psum(i),i=1,maxp) 1164 write(lfnmro) (p2sum(i),i=1,maxp) 1165 write(lfnmro) (pslop(i),i=1,maxp) 1166 write(lfnmro) (psumt(i),i=1,maxp) 1167 write(lfnmro) (p2sumt(i),i=1,maxp) 1168 write(lfnmro) (pslopt(i),i=1,maxp) 1169 endif 1170 if(ndec.gt.0) call cf_wrtmro(lfnmro) 1171c 1172 return 1173 end 1174 logical function prp_rmri(lfnmri,ndec,mropt,dfr,dfrm) 1175c 1176 implicit none 1177c 1178#include "prp_common.fh" 1179c 1180 logical cf_rdmri 1181 external cf_rdmri 1182c 1183 integer lfnmri,ndec,mropt 1184 real*8 dfr(mda),dfrm(mda),dtmp 1185c 1186 integer i,nprp,nmxp 1187c 1188 read(lfnmri,err=9,end=9) nderiv,nda,nprp,nsum,nsumt,ndec,nmxp 1189 if(mropt.ne.2) then 1190 if(nprop.ne.nprp) then 1191 call md_abort('Number of properties nprop changed',0) 1192 endif 1193 if(maxp.ne.nmxp) then 1194 call md_abort('Number of properties maxp changed',0) 1195 endif 1196 endif 1197 read(lfnmri,err=9,end=9) deriv 1198 if(mda.lt.nda) then 1199 read(lfnmri,err=9,end=9) (dtmp,i=1,nda) 1200 read(lfnmri,err=9,end=9) (dtmp,i=1,nda) 1201 else 1202 read(lfnmri,err=9,end=9) (dfr(i),i=1,nda) 1203 read(lfnmri,err=9,end=9) (dfrm(i),i=1,nda) 1204 endif 1205 read(lfnmri,err=9,end=9) tsum,t2sum,tsumt,t2sumt 1206 if(mropt.ne.2) then 1207 read(lfnmri,err=9,end=9) (psum(i),i=1,maxp) 1208 read(lfnmri,err=9,end=9) (p2sum(i),i=1,maxp) 1209 read(lfnmri,err=9,end=9) (pslop(i),i=1,maxp) 1210 read(lfnmri,err=9,end=9) (psumt(i),i=1,maxp) 1211 read(lfnmri,err=9,end=9) (p2sumt(i),i=1,maxp) 1212 read(lfnmri,err=9,end=9) (pslopt(i),i=1,maxp) 1213 else 1214 do 1 i=1,6 1215 read(lfnmri,err=9,end=9) dtmp 1216 1 continue 1217 endif 1218 if(ndec.gt.0) then 1219 prp_rmri=cf_rdmri(lfnmri) 1220 else 1221 prp_rmri=.true. 1222 endif 1223 return 1224c 1225 9 continue 1226 prp_rmri=.false. 1227 return 1228 end 1229 integer function prp_dfr(ncopy) 1230c 1231 implicit none 1232c 1233#include "prp_common.fh" 1234#include "mafdecls.fh" 1235c 1236 integer prp_dfr_copy 1237 external prp_dfr_copy 1238c 1239 integer ncopy 1240c 1241 prp_dfr=prp_dfr_copy(ncopy,dbl_mb(i_dfr),dbl_mb(i_dfrm)) 1242c 1243 return 1244 end 1245 integer function prp_dfr_copy(ncopy,dfr,dfrm) 1246c 1247 implicit none 1248c 1249#include "prp_common.fh" 1250c 1251 integer ncopy 1252 real*8 dfr(mda),dfrm(mda) 1253c 1254 integer i,ioff 1255c 1256 if(ncopy.lt.nda) then 1257 ioff=nda-ncopy 1258 do 1 i=1,ncopy 1259 dfr(i)=dfr(i+ioff) 1260 dfrm(i)=dfrm(i+ioff) 1261 1 continue 1262 nda=ncopy 1263 endif 1264c 1265 prp_dfr_copy=nda 1266 return 1267 end 1268