1 subroutine md_main() 2c 3c $Id$ 4c 5 implicit none 6c 7#include "md_common.fh" 8c 9 if(ntype.eq.0) then 10c 11c single energy 12c ------------- 13c 14 if(nftri.eq.0) then 15 call md_sp() 16 else 17 call md_spi() 18 endif 19c 20 elseif(ntype.eq.1) then 21c 22c energy minimization 23c ------------------- 24c 25 call md_em() 26c 27 elseif(ntype.eq.2) then 28c 29c molecular dynamics 30c ------------------ 31c 32 call md_md() 33c 34 elseif(ntype.eq.3) then 35c 36c free energy simulation 37c ---------------------- 38c 39 call md_ti() 40c 41 else 42 call md_abort('Unknown calculation type',ntype) 43 endif 44c 45 return 46 end 47 subroutine md_sp() 48c 49 implicit none 50c 51#include "md_common.fh" 52#include "mafdecls.fh" 53c 54 lpair=.true. 55 lload=.false. 56 lhop=.false. 57 llong=ltwin 58c 59c center of mass coordinates 60c 61 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 62 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 63 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 64c 65c periodic boundary conditions 66c 67 call md_fold(int_mb(i_iw),int_mb(i_is), 68 + dbl_mb(i_xw),dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_xsm)) 69c 70c atom redistribution 71c 72 call sp_travel(box,dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 73 + dbl_mb(i_gw),int_mb(i_iw),nwmloc,dbl_mb(i_xs),dbl_mb(i_vs), 74 + dbl_mb(i_gs),int_mb(i_is),nsaloc) 75c 76c center of mass coordinates 77c 78 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 79 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 80 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 81c 82 call cf_mass(dbl_mb(i_wws),dbl_mb(i_wws+mwa), 83 + int_mb(i_is+(lsatt-1)*msa),nsaloc) 84c 85 call md_eminit(dbl_mb(i_xw),dbl_mb(i_yw), 86 + dbl_mb(i_xs),dbl_mb(i_ys)) 87c 88c atomic forces and potential energies 89c 90 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 91 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 92 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 93 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 94 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 95c 96 call prp_proper(0,stime,eww,dbl_mb(i_esw), 97 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 98 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 99 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 100c 101 if(.not.lqmmm) call prp_print() 102c 103 call rtdb_put(irtdb,'md:energy',mt_dbl,1,epot) 104c 105c print energies 106c 107 if(.not.lqmmm) then 108 call cf_print_energy(lfnout) 109 call sp_printf(filtop,lfntop, 110 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_fs),npener,dbl_mb(i_esa)) 111 endif 112c 113 if(ifidi.ne.0) then 114 call md_fd(int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),dbl_mb(i_fs), 115 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw)) 116 endif 117c 118 if(itest.eq.1) call md_test() 119c 120 return 121 end 122 subroutine md_sp_qmmm() 123c 124 implicit none 125c 126#include "md_common.fh" 127#include "mafdecls.fh" 128#include "msgids.fh" 129 lpair=.false. 130 lload=.false. 131 132c atomic forces and potential energies 133c 134 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 135 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 136 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 137 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 138 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 139c 140 call prp_proper(0,stime,eww,dbl_mb(i_esw), 141 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 142 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 143 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 144 145 call rtdb_put(irtdb,'md:energy',mt_dbl,1,epot) 146 147 return 148 end 149 subroutine md_spi() 150c 151 implicit none 152c 153#include "md_common.fh" 154#include "mafdecls.fh" 155c 156 external sp_rdtrj,sp_skip,frequency 157 logical sp_rdtrj,sp_skip,frequency 158c 159 integer import 160c 161 if(impfr.gt.1) then 162 if(.not.sp_skip(lfntri,impfr-1)) return 163 endif 164c 165 do 1 import=impfr,impto 166c 167 if(.not.frequency(import+1-impfr,nftri)) then 168 if(.not.sp_skip(lfntri,1)) return 169 goto 1 170 endif 171c 172 lpair=.true. 173 lload=.true. 174 lhop=.false. 175 llong=ltwin 176c 177c read frame from trajectory file 178c 179 if(.not.sp_rdtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs, 180 + stime,pres,temp,tempw,temps, 181 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw), 182 + dbl_mb(i_xwcr),int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs), 183 + dbl_mb(i_fs),nwmloc,nsaloc)) return 184c 185c center of mass coordinates 186c 187 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 188 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 189 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 190c 191c periodic boundary conditions 192c 193 call md_fold(int_mb(i_iw),int_mb(i_is), 194 + dbl_mb(i_xw),dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_xsm)) 195c 196c atom redistribution 197c 198 call sp_travel(box,dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 199 + dbl_mb(i_gw),int_mb(i_iw),nwmloc,dbl_mb(i_xs),dbl_mb(i_vs), 200 + dbl_mb(i_gs),int_mb(i_is),nsaloc) 201c 202 call cf_mass(dbl_mb(i_wws),dbl_mb(i_wws+mwa), 203 + int_mb(i_is+(lsatt-1)*msa),nsaloc) 204c 205 call md_eminit(dbl_mb(i_xw),dbl_mb(i_yw), 206 + dbl_mb(i_xs),dbl_mb(i_ys)) 207c 208c atomic forces and potential energies 209c 210 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 211 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 212 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 213 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 214 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 215c 216 call prp_proper(import,stime,eww,dbl_mb(i_esw), 217 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 218 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 219 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 220c 221 call prp_step(import,stime,eww,dbl_mb(i_esw), 222 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm) 223c 224cx call prp_print() 225c 226 call rtdb_put(irtdb,'md:energy',mt_dbl,1,epot) 227c 228c print energies 229c 230 if(me.eq.0.and.frequency(import,npener)) 231 + call cf_print_energy(lfnout) 232c if(me.eq.0.and.frequency(import,nfprop)) call prp_record() 233c 234 if(me.eq.0.and.frequency(import,npforc)) 235 + call sp_printf(filtop,lfntop, 236 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_fs),npener,dbl_mb(i_esa)) 237c 238 if(ifidi.ne.0) then 239 call md_fd(int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),dbl_mb(i_fs), 240 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw)) 241 endif 242c 243 if(itest.eq.1) call md_test() 244c 245 1 continue 246c 247 return 248 end 249 subroutine md_em() 250c 251 implicit none 252c 253#include "md_common.fh" 254#include "mafdecls.fh" 255c 256 integer i_pcgw,l_pcgw,i_pcgs,l_pcgs 257c 258 llong=ltwin 259c 260c center of mass coordinates 261c 262 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 263 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 264 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 265c 266c periodic boundary conditions 267c 268 call md_fold(int_mb(i_iw),int_mb(i_is), 269 + dbl_mb(i_xw),dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_xsm)) 270c 271c atom redistribution 272c 273 call sp_travel(box,dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 274 + dbl_mb(i_gw),int_mb(i_iw),nwmloc,dbl_mb(i_xs),dbl_mb(i_vs), 275 + dbl_mb(i_gs),int_mb(i_is),nsaloc) 276c 277c center of mass coordinates 278c 279 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 280 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 281 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 282c 283 call cf_mass(dbl_mb(i_wws),dbl_mb(i_wws+mwa), 284 + int_mb(i_is+(lsatt-1)*msa),nsaloc) 285c 286 call md_eminit(dbl_mb(i_xw),dbl_mb(i_yw), 287 + dbl_mb(i_xs),dbl_mb(i_ys)) 288c 289 if(msdit.gt.0) call md_stdesc(int_mb(i_iw+(lwdyn-1)*mwm), 290 + dbl_mb(i_xw),dbl_mb(i_yw),dbl_mb(i_vw),dbl_mb(i_fw), 291 + int_mb(i_is+(lsdyn-1)*msa),int_mb(i_is+(lsmol-1)*msa), 292 + dbl_mb(i_xs),dbl_mb(i_ys),dbl_mb(i_vs),dbl_mb(i_fs), 293 + dbl_mb(i_wws),dbl_mb(i_wws+mwa),int_mb(i_mm),dbl_mb(i_fm), 294 + dbl_mb(i_xsm)) 295c 296 if(mcgit.gt.0) then 297 if(.not.ma_push_get(mt_dbl,3*mwa*mwm,'pcgw',l_pcgw,i_pcgw)) 298 + call md_abort('Failed to allocate memory for pcgw',me) 299 if(.not.ma_push_get(mt_dbl,3*msa,'pcgs',l_pcgs,i_pcgs)) 300 + call md_abort('Failed to allocate memory for pcgs',me) 301 call md_congra(int_mb(i_iw+(lwdyn-1)*mwm),dbl_mb(i_xw), 302 + dbl_mb(i_yw),dbl_mb(i_vw),dbl_mb(i_fw),dbl_mb(i_pcgw), 303 + int_mb(i_is+(lsdyn-1)*msa),dbl_mb(i_xs),dbl_mb(i_ys), 304 + dbl_mb(i_vs),dbl_mb(i_fs),dbl_mb(i_pcgs),dbl_mb(i_wws), 305 + dbl_mb(i_wws+mwa)) 306 if(.not.ma_pop_stack(l_pcgs)) 307 + call md_abort('Failed to deallocate memory for pcgs',me) 308 if(.not.ma_pop_stack(l_pcgw)) 309 + call md_abort('Failed to deallocate memory for pcgw',me) 310 endif 311c 312 call cf_print_energy(lfnout) 313c 314 call sp_printf(filtop,lfntop, 315 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_fs),npener,dbl_mb(i_esa)) 316c 317 return 318 end 319 subroutine md_eminit(xw,yw,xs,ys) 320c 321 implicit none 322c 323#include "md_common.fh" 324#include "mafdecls.fh" 325c 326 real*8 xw(mwm,3,mwa),yw(mwm,3,mwa) 327 real*8 xs(msa,3),ys(msa,3) 328c 329 integer i,j 330 real*8 dxmax 331c 332 if(nwmloc.gt.0) then 333 do 1 j=1,nwa 334 do 2 i=1,nwmloc 335 yw(i,1,j)=xw(i,1,j) 336 yw(i,2,j)=xw(i,2,j) 337 yw(i,3,j)=xw(i,3,j) 338 2 continue 339 1 continue 340 endif 341c 342 if(nsaloc.gt.0) then 343 do 3 i=1,nsaloc 344 ys(i,1)=xs(i,1) 345 ys(i,2)=xs(i,2) 346 ys(i,3)=xs(i,3) 347 3 continue 348 endif 349c 350 return 351c 352 lpair=.true. 353 lload=.true. 354 lhop=.false. 355c 356 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 357 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 358 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 359 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 360 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 361c 362c shake 363c 364 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 365 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 366c 367 return 368 end 369 subroutine md_stdesc(iwdt,xw,yw,vw,fw,isdt,ismol, 370 + xs,ys,vs,fs,ww,ws,mm,fm,xsm) 371c 372 implicit none 373c 374#include "md_common.fh" 375#include "mafdecls.fh" 376#include "msgids.fh" 377#include "global.fh" 378c 379 logical frequency 380 external frequency 381c 382 integer iwdt(mwm),isdt(msa),ismol(msa),mm(msa,2) 383 real*8 xw(mwm,3,mwa),yw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa) 384 real*8 xs(msa,3),ys(msa,3),vs(msa,3),fs(msa,3) 385 real*8 ww(mwa),ws(msa),xsm(msm,3),fm(msm,7) 386c 387 integer i,j 388 logical ldone 389 real*8 edif,epsd,epsdw,epsdsw,epsds,ecrit 390 real*8 da,damsd,dx,dxf,dxmax,factor,dxstep,eqrs 391 character*1 cqrs 392 real*8 angle,o(3),p(3),x(3),y(3) 393c 394 double precision grms(2) 395 integer nrms(2) 396c 397 damsd=-1.1d0 398 if(imembr.gt.0) call md_membrane_init(ismol,mm,xs,xsm,fm) 399c 400 if(lqmmm) then 401 if(me.eq.0) write(6,601) "@","@" 402 end if 403601 format( 404 $ /,a1,' Step Energy Delta E SGrms', 405 $ ' WGrms Xmax', 406 $ /,a1,' ---- ---------------- -------- --------', 407 $ ' -------- -------- ') 408 409 if(me.eq.0) write(lfnout,1000) 410 1000 format(/,' STEEPEST DESCENT MINIMIZATION',//, 411 + ' Step File Energy Energy Energy ', 412 + ' Energy Energy Largest ',/, 413 + ' wrt gradient Total solvent ', 414 + ' slv-sol solute displacement',/, 415 + ' kJ/mol kJ/mol kJ/mol ', 416 + ' kJ/mol kJ/mol nm',/) 417c 418 isdit=0 419c 420 lpair=.true. 421 lload=.true. 422 lhop=.false. 423c 424c atomic forces and potential energies 425c 426 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 427 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 428 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 429 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 430 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 431c 432 call prp_proper(isdit,stime,eww,dbl_mb(i_esw), 433 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 434 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 435 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 436c 437 if(lqmmm) call qmmm_print_energy1(irtdb,epot,uqmmm) 438c 439 epsd=epot 440 epsdw=epotw 441 epsdsw=epotsw 442 epsds=epots 443 eqrs=epot 444 dx=dx0sd 445 mdstep=0 446 da=0.01 447c 448 1 continue 449 call timer_start(201) 450c 451 dxf=dx/fmax 452c 453 if(nwmloc.gt.0) then 454 do 2 j=1,nwa 455 do 3 i=1,nwmloc 456 yw(i,1,j)=xw(i,1,j) 457 yw(i,2,j)=xw(i,2,j) 458 yw(i,3,j)=xw(i,3,j) 459 vw(i,1,j)=fw(i,1,j) 460 vw(i,2,j)=fw(i,2,j) 461 vw(i,3,j)=fw(i,3,j) 462 3 continue 463 factor=one/ww(j) 464 do 4 i=1,nwmloc 465 if(iand(iwdt(i),mfixed).ne.lfixed) then 466 dxstep=factor*dxf*fw(i,1,j) 467 xw(i,1,j)=xw(i,1,j)+dxstep 468 dxstep=factor*dxf*fw(i,2,j) 469 xw(i,2,j)=xw(i,2,j)+dxstep 470 dxstep=factor*dxf*fw(i,3,j) 471 xw(i,3,j)=xw(i,3,j)+dxstep 472 endif 473 4 continue 474 2 continue 475 endif 476c 477 if(nsaloc.gt.0) then 478 do 5 i=1,nsaloc 479 ys(i,1)=xs(i,1) 480 ys(i,2)=xs(i,2) 481 ys(i,3)=xs(i,3) 482 vs(i,1)=fs(i,1) 483 vs(i,2)=fs(i,2) 484 vs(i,3)=fs(i,3) 485 5 continue 486 if(imembr.eq.0) then 487 do 6 i=1,nsaloc 488 if(iand(isdt(i),mfixed).ne.lfixed) then 489 factor=one/ws(i) 490 dxstep=factor*dxf*fs(i,1) 491 xs(i,1)=xs(i,1)+dxstep 492 dxstep=factor*dxf*fs(i,2) 493 xs(i,2)=xs(i,2)+dxstep 494 dxstep=factor*dxf*fs(i,3) 495 xs(i,3)=xs(i,3)+dxstep 496 endif 497 6 continue 498 else 499 do 16 i=1,msm 500 fm(i,1)=zero 501 fm(i,2)=zero 502 fm(i,3)=zero 503 fm(i,4)=zero 504 fm(i,5)=zero 505 fm(i,6)=zero 506 fm(i,7)=zero 507 16 continue 508 do 17 i=1,nsaloc 509 factor=one/ws(i) 510 fm(mm(i,2),1)=fm(mm(i,2),1)+factor*fs(i,1) 511 fm(mm(i,2),2)=fm(mm(i,2),2)+factor*fs(i,2) 512 fm(mm(i,2),3)=fm(mm(i,2),3)+factor* 513 + ((xs(i,1)-xsm(mm(i,2),1))*fs(i,2)- 514 + (xs(i,2)-xsm(mm(i,2),2))*fs(i,1)) 515 17 continue 516 if(np.gt.1) call ga_dgop(mrg_d50,fm,3*msm,'+') 517 do 18 i=1,nsm 518 fm(i,4)=fm(i,3) 519 write(*,'(a,i5,a,i5)') 'Mol ',i,', Angle ',da*fm(i,4) 520 18 continue 521 do 19 i=1,nsaloc 522 if(iand(isdt(i),mfixed).ne.lfixed) then 523c 524c rotations 525c 526 o(1)=xsm(mm(i,2),1) 527 o(2)=xsm(mm(i,2),1) 528 o(3)=xsm(mm(i,2),1) 529 p(1)=o(1) 530 p(2)=o(2) 531 p(3)=o(3)+1.0d0 532 x(1)=xs(i,1) 533 x(2)=xs(i,2) 534 x(3)=xs(i,3) 535 angle=da*fm(mm(i,2),4) 536 call rotate(o,p,angle,x,y) 537 xs(i,1)=y(1) 538 xs(i,2)=y(2) 539 xs(i,3)=y(3) 540c 541c translations 542c 543 dxstep=dxf*fm(mm(i,2),1) 544 xs(i,1)=xs(i,1)+dxstep 545 dxstep=dxf*fm(mm(i,2),2) 546 xs(i,2)=xs(i,2)+dxstep 547c 548 endif 549 19 continue 550 endif 551 endif 552c 553c shake 554c 555 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 556 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 557c 558 isdit=isdit+1 559c 560 lpair=frequency(isdit,nfpair) 561 lload=lpair 562 lhop=.false. 563 llong=(frequency(isdit,nflong).or.lpair).and.ltwin 564c 565c atomic forces and potential energies 566c 567 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 568 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 569 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 570 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 571 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 572c 573 call prp_proper(isdit,stime,eww,dbl_mb(i_esw), 574 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 575 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 576 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 577c 578 if(lqmmm) call qmmm_print_energy1(irtdb,epot,uqmmm) 579c 580c if energy goes up restore coordinates and forces 581c 582 ecrit=epot 583 if(icrit.eq.1) ecrit=epot 584c 585c changing the logic here so that bad coordinates 586c are not preserved in the last md iteration (M.V.) 587c ------------------------------------------------ 588 if(ecrit.gt.epsd) then 589c 590 call cf_restor(dbl_mb(i_xw),dbl_mb(i_yw),dbl_mb(i_fw), 591 + dbl_mb(i_vw),nwmloc,dbl_mb(i_xs),dbl_mb(i_ys),dbl_mb(i_fs), 592 + dbl_mb(i_vs),nsaloc) 593c 594 if(isdit.lt.msdit.and.dxmax.gt.dxsdmx) then 595 dx=half*dx 596 da=half*da 597 call timer_stop(201) 598 if(dx.gt.dxsdmx) goto 1 599 else 600c 601c this insures the global restoration of coordinates 602c so that the right restart file written out 603c -------------------------------------------------- 604 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 605 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 606 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 607 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 608 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 609c 610 call prp_proper(isdit,stime,eww,dbl_mb(i_esw), 611 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 612 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 613 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 614c 615 endif 616 endif 617c 618 edif=ecrit-epsd 619 epsd=ecrit 620 epsdw=epotw 621 epsdsw=epotsw 622 epsds=epots 623c 624 lxw=frequency(mdstep,nfcoor) 625 lxs=frequency(mdstep,nfscoo) 626c 627 if(lxw.or.lxs) then 628 if(lqmd) then 629 call qmd_wrttrj(lfntrj,mwm,nwmloc,mwa,nwa,msa,nsaloc, 630 + .true.,.false.,.true.,.false.,box,stime,pres,temp,tempw,temps, 631 + dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xs),dbl_mb(i_vs)) 632 else 633 call sp_wrttrj(lfntrj,lxw,.false.,.false.,lxs,.false.,.false., 634 + stime,pres,temp,tempw,temps, 635 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw), 636 + dbl_mb(i_xwcr),int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs), 637 + dbl_mb(i_fs)) 638 endif 639 endif 640c 641 ldone=.not.(dxmax.gt.dxsdmx.and.isdit.lt.msdit) 642c 643 cqrs=' ' 644 if(frequency(isdit,nfqrs).or.(ldone.and.eqrs.gt.epot)) then 645 write(projct,4000) nserie,isdit,0,filnam(1:56) 646 4000 format(i2,' em ',i7,' + ',i7,' ',a) 647 call md_wrtrst(lfnqrs,filqrs,.false.) 648 cqrs='*' 649 eqrs=epot 650 endif 651 if(me.eq.0.and.frequency(isdit,nfprop)) call prp_record() 652c 653 if(me.eq.0) then 654 write(lfnout,600) isdit,cqrs,edif,epsd,epsdw,epsdsw,epsds,dxmax 655 600 format(i7,1x,a1,3x,5(1pe13.5),0pf12.8) 656 call util_flush(lfnout) 657 endif 658c 659 if(lqmmm) then 660 nrms(2) = 0 661 grms(2) = 0.0d0 662 if(nwmloc.gt.0) then 663 do j=1,nwa 664 do i=1,nwmloc 665 if(iand(iwdt(i),mfixed).ne.lfixed) then 666 nrms(2) = nrms(2) + 1 667 grms(2) = grms(2) + 668 + fw(i,1,j)*fw(i,1,j) + 669 + fw(i,2,j)*fw(i,2,j) + 670 + fw(i,3,j)*fw(i,3,j) 671 endif 672 end do 673 end do 674 endif 675c 676 nrms(1) = 0 677 grms(1) = 0.0d0 678 do i=1,nsaloc 679 if(iand(isdt(i),mfixed).ne.lfixed) then 680 nrms(1) = nrms(1) + 1 681 grms(1) = grms(1) + 682 + fs(i,1)*fs(i,1)+ 683 + fs(i,2)*fs(i,2)+ 684 + fs(i,3)*fs(i,3) 685 endif 686 end do 687 call ga_dgop(msg_qmmm_force,grms,2,'+') 688 call ga_igop(msg_qmmm_nat1,nrms,2,'+') 689 690 do i=1,2 691 if(nrms(i).ne.0) then 692 grms(i) = sqrt(grms(i)/dble(nrms(i))) 693 grms(i) = grms(i)*5.29177249d-02/2.625499962d+03 694 end if 695 end do 696 if(me.eq.0) 697 $ write(6,602) "@", isdit, epsd/2.625499962d+03, 698 $ edif/2.625499962d+03, 699 $ grms(1), grms(2), 700 $ dxmax/5.29177249d-02 701602 format(a1,i5,f17.8,1p,d9.1,0p,3f9.5) 702 703 end if 704c 705 706 if(itest.gt.0) call md_test() 707c 708 dx=min(1.2d0*dx,dxmsd) 709 da=min(1.2d0*da,damsd) 710c 711 call timer_stop(201) 712 if(.not.ldone) goto 1 713 714 return 715 end 716 subroutine md_congra(iwdt,xw,yw,vw,fw,pcgw, 717 + isdt,xs,ys,vs,fs,pcgs,ww,ws) 718c 719 implicit none 720c 721#include "md_common.fh" 722#include "mafdecls.fh" 723#include "msgids.fh" 724c 725 logical frequency 726 external frequency 727c 728 integer iwdt(mwm),isdt(msa) 729 real*8 xw(mwm,3,mwa),yw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa) 730 real*8 xs(msa,3),ys(msa,3),vs(msa,3),fs(msa,3) 731 real*8 pcgw(mwm,3,mwa),pcgs(msa,3),ww(mwa),ws(msa) 732c 733 integer iwa,iwm,isa,ix,inner 734 logical ldone 735 real*8 alpha,beta1,beta2,beta3,beta4,beta5,gamma,zeta 736 real*8 dx,ecgnew,eqrs,ecgold,ecgdif,edt,dxf 737 real*8 dxsmax,dxwmax,fnorm1,fnorm2,ypa,ypb,pnorm,dxstep 738 real*8 dxmax,dxfi,ecg0,ecg1,ecg2,epcgw,epcgsw,epcgs 739 character*1 cqrs 740c 741 if(me.eq.0) write(lfnout,1000) 742 1000 format(/,' CONJUGATE GRADIENT MINIMIZATION',//, 743 + ' Step File Energy Energy Energy ', 744 + ' Energy Energy Largest ',/, 745 + ' wrt gradient Total solvent ', 746 + ' slv-sol solute displacement',/, 747 + ' kJ/mol kJ/mol kJ/mol ', 748 + ' kJ/mol kJ/mol nm',/) 749c 750 dx=dx0cg 751 beta1=zero 752 icgit=0 753 lpair=.true. 754 lload=.true. 755 lhop=.false. 756c 757 call timer_start(201) 758c 759c atomic forces and potential energies 760c 761 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 762 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 763 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 764 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 765 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 766c 767 call prp_proper(isdit,stime,eww,dbl_mb(i_esw), 768 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 769 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 770 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 771c 772 ecgnew=epot 773 eqrs=epot 774 ecgold=ecgnew 775 ecgdif=zero 776 edt=zero 777c 778c copy initial coordinates into vw and vs 779c initialize search direction vectors pcgw and pcgs to zero 780c 781 if(nwmloc.gt.0) then 782 do 1 iwa=1,nwa 783 do 2 ix=1,3 784 do 3 iwm=1,nwmloc 785 vw(iwm,ix,iwa)=xw(iwm,ix,iwa) 786 pcgw(iwm,ix,iwa)=zero 787 3 continue 788 2 continue 789 1 continue 790 endif 791 if(nsaloc.gt.0) then 792 do 4 ix=1,3 793 do 5 isa=1,nsaloc 794 vs(isa,ix)=xs(isa,ix) 795 pcgs(isa,ix)=zero 796 5 continue 797 4 continue 798 endif 799c 800c take one steepest descent step to get initial search direction 801c 802 dxf=half*dx/fmax 803 if(nwmloc.gt.0) then 804 do 6 iwa=1,nwa 805 do 7 ix=1,3 806 do 8 iwm=1,nwmloc 807 yw(iwm,ix,iwa)=xw(iwm,ix,iwa) 808 if(iand(iwdt(iwm),mfixed).ne.lfixed) then 809 xw(iwm,ix,iwa)=xw(iwm,ix,iwa)+dxf*fw(iwm,ix,iwa)/ww(iwa) 810 endif 811 8 continue 812 7 continue 813 6 continue 814 endif 815 if(nsaloc.gt.0) then 816 do 9 ix=1,3 817 do 10 isa=1,nsaloc 818 ys(isa,ix)=xs(isa,ix) 819 if(iand(isdt(isa),mfixed).ne.lfixed) then 820 xs(isa,ix)=xs(isa,ix)+dxf*fs(isa,ix)/ws(isa) 821 endif 822 10 continue 823 9 continue 824 endif 825c 826c shake 827c 828 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 829 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 830c 831 fnorm1=zero 832 dxfi=one/dxf 833 if(nwmloc.gt.0) then 834 do 11 iwa=1,nwa 835 do 12 ix=1,3 836 do 13 iwm=1,nwmloc 837 fw(iwm,ix,iwa)=(xw(iwm,ix,iwa)-yw(iwm,ix,iwa))*dxfi*ww(iwa) 838 fnorm1=fnorm1+fw(iwm,ix,iwa)**2 839 13 continue 840 12 continue 841 11 continue 842 endif 843 if(nsaloc.gt.0) then 844 do 14 ix=1,3 845 do 15 isa=1,nsaloc 846 fs(isa,ix)=(xs(isa,ix)-ys(isa,ix))*dxfi*ws(isa) 847 fnorm1=fnorm1+fs(isa,ix)**2 848 15 continue 849 14 continue 850 endif 851c 852c global sum fnorm1 853c 854 call ga_dgop(mrg_d08,fnorm1,1,'+') 855c 856 ecg0=ecgnew 857 ecg1=ecgnew 858 ecg2=ecgnew 859 icgit=0 860c 861 call timer_stop(201) 862c 863c outer loop 864c 865 100 continue 866c 867 if(icgit.eq.(icgit/ncgcy)*ncgcy) beta1=zero 868 icgit=icgit+1 869 lpair=frequency(icgit,nfpair) 870 lload=frequency(icgit,nfload) 871 lhop=.false. 872c 873 ypa=zero 874 pnorm=zero 875 if(nwmloc.gt.0) then 876 do 16 iwa=1,nwa 877 do 17 ix=1,3 878 do 18 iwm=1,nwmloc 879 pcgw(iwm,ix,iwa)=fw(iwm,ix,iwa)+beta1*pcgw(iwm,ix,iwa) 880 pnorm=pnorm+pcgw(iwm,ix,iwa)**2 881 ypa=ypa+pcgw(iwm,ix,iwa)*fw(iwm,ix,iwa) 882 18 continue 883 17 continue 884 16 continue 885 endif 886 if(nsaloc.gt.0) then 887 do 19 ix=1,3 888 do 20 isa=1,nsaloc 889 pcgs(isa,ix)=fs(isa,ix)+beta1*pcgs(isa,ix) 890 pnorm=pnorm+pcgs(isa,ix)**2 891 ypa=ypa+pcgs(isa,ix)*fs(isa,ix) 892 20 continue 893 19 continue 894 endif 895c 896c accumulate pnorm 897c 898 call ga_dgop(mrg_d09,ypa,1,'+') 899 call ga_dgop(mrg_d10,pnorm,1,'+') 900c 901 if(pnorm.lt.zero) call md_abort('congra: pnorm<zero',0) 902 pnorm=sqrt(pnorm) 903c 904 alpha=zero 905 ecg1=ecg2 906 beta2=dx/pnorm 907 inner=0 908c 909c inner loop 910c 911 200 continue 912 call timer_start(201) 913c 914 inner=inner+1 915c 916 if(nwmloc.gt.0) then 917 do 21 iwa=1,nwa 918 do 22 ix=1,3 919 do 23 iwm=1,nwmloc 920 if(iand(iwdt(iwm),mfixed).ne.lfixed) then 921 xw(iwm,ix,iwa)=vw(iwm,ix,iwa)+beta2*pcgw(iwm,ix,iwa)/ww(iwa) 922 endif 923 23 continue 924 22 continue 925 21 continue 926 endif 927 if(nsaloc.gt.0) then 928 do 24 ix=1,3 929 do 25 isa=1,nsaloc 930 if(iand(isdt(isa),mfixed).ne.lfixed) then 931 xs(isa,ix)=vs(isa,ix)+beta2*pcgs(isa,ix)/ws(isa) 932 endif 933 25 continue 934 24 continue 935 endif 936c 937c atomic forces and potential energies 938c 939 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 940 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 941 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 942 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 943 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 944c 945 lpair=.false. 946 lload=.false. 947 lhop=.false. 948c 949 beta3=beta2*pnorm/fnorm 950c 951 if(nwmloc.gt.0) then 952 do 26 iwa=1,nwa 953 do 27 ix=1,3 954 do 28 iwm=1,nwmloc 955 yw(iwm,ix,iwa)=xw(iwm,ix,iwa) 956 dxstep=beta3*fw(iwm,ix,iwa)/ww(iwa) 957 if(iand(iwdt(iwm),mfixed).ne.lfixed) then 958 xw(iwm,ix,iwa)=yw(iwm,ix,iwa)+dxstep 959 endif 960 28 continue 961 27 continue 962 26 continue 963 endif 964 if(nsaloc.gt.0) then 965 do 29 ix=1,3 966 do 30 isa=1,nsaloc 967 ys(isa,ix)=xs(isa,ix) 968 dxstep=beta3*fs(isa,ix)/ws(isa) 969 if(iand(isdt(isa),mfixed).ne.lfixed) then 970 xs(isa,ix)=ys(isa,ix)+dxstep 971 endif 972 30 continue 973 29 continue 974 endif 975c 976c shake 977c 978 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 979 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 980c 981c find constrained forces 982c 983 ypb=zero 984 if(nwmloc.gt.0) then 985 do 34 iwa=1,nwa 986 do 35 ix=1,3 987 do 36 iwm=1,nwmloc 988 fw(iwm,ix,iwa)=(xw(iwm,ix,iwa)-yw(iwm,ix,iwa))*ww(iwa)/beta3 989 ypb=ypb+pcgw(iwm,ix,iwa)*fw(iwm,ix,iwa) 990 36 continue 991 35 continue 992 34 continue 993 endif 994 if(nsaloc.gt.0) then 995 do 37 ix=1,3 996 do 38 isa=1,nsaloc 997 fs(isa,ix)=(xs(isa,ix)-ys(isa,ix))*ws(isa)/beta3 998 ypb=ypb+pcgs(isa,ix)*fs(isa,ix) 999 38 continue 1000 37 continue 1001 endif 1002 call ga_dgop(mrg_d11,ypb,1,'+') 1003c 1004 call prp_proper(isdit,stime,eww,dbl_mb(i_esw), 1005 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 1006 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 1007 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 1008c 1009 ecg2=epot 1010 epcgw=epotw 1011 epcgsw=epotsw 1012 epcgs=epots 1013 ecgdif=ecg2-ecgold 1014 edt=ecg2-ecg0 1015 call timer_stop(201) 1016c 1017c check if interval is appropriate 1018c 1019 if(ypb.ge.zero.and.ecg2.lt.ecg1) then 1020 alpha=beta2 1021 ecg1=ecg2 1022 ypa=ypb 1023 beta2=two*beta2 1024 goto 200 1025 endif 1026 call timer_start(201) 1027c 1028c find interpolation in interval 1029c 1030 zeta=three*(ecg1-ecg2)/(beta2-alpha)-ypa-ypb 1031 gamma=zeta**2-ypa*ypb 1032c 1033 if(gamma.lt.zero) then 1034 gamma=zero 1035 else 1036 gamma=sqrt(gamma) 1037 endif 1038c 1039 beta4=beta2-(gamma-zeta-ypb)*(beta2-alpha)/(ypa-ypb+two*gamma) 1040c 1041c advance coordinates to interpolated point 1042c 1043 dxmax=zero 1044 if(nwmloc.gt.0) then 1045 do 39 iwa=1,nwa 1046 do 40 ix=1,3 1047 do 41 iwm=1,nwmloc 1048 yw(iwm,ix,iwa)=vw(iwm,ix,iwa) 1049 dxstep=beta4*pcgw(iwm,ix,iwa)/ww(iwa) 1050 if(iand(iwdt(iwm),mfixed).ne.lfixed) then 1051 xw(iwm,ix,iwa)=vw(iwm,ix,iwa)+dxstep 1052 if(abs(dxstep).gt.dxmax) dxmax=abs(dxstep) 1053 endif 1054 41 continue 1055 40 continue 1056 39 continue 1057 endif 1058 if(nsaloc.gt.0) then 1059 do 42 ix=1,3 1060 do 43 isa=1,nsaloc 1061 ys(isa,ix)=vs(isa,ix) 1062 dxstep=beta4*pcgs(isa,ix)/ws(isa) 1063 if(iand(isdt(isa),mfixed).ne.lfixed) then 1064 xs(isa,ix)=vs(isa,ix)+dxstep 1065 if(abs(dxstep).gt.dxmax) dxmax=abs(dxstep) 1066 endif 1067 43 continue 1068 42 continue 1069 endif 1070 call ga_dgop(mrg_d12,dxmax,1,'max') 1071c 1072c shake 1073c 1074 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 1075 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 1076c 1077c atomic forces and potential energies 1078c 1079 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 1080 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 1081 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 1082 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 1083 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 1084c 1085 call prp_proper(isdit,stime,eww,dbl_mb(i_esw), 1086 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 1087 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 1088 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 1089c 1090 dxmax=zero 1091 dxwmax=zero 1092 dxsmax=zero 1093c 1094 beta5=beta4*pnorm/fnorm 1095c 1096c advance coordinates with these forces 1097c 1098 if(nwmloc.gt.0) then 1099 do 44 iwa=1,nwa 1100 do 45 ix=1,3 1101 do 46 iwm=1,nwmloc 1102 yw(iwm,ix,iwa)=xw(iwm,ix,iwa) 1103 dxstep=beta5*fw(iwm,ix,iwa)/ww(iwa) 1104 if(iand(iwdt(iwm),mfixed).ne.lfixed) then 1105 xw(iwm,ix,iwa)=yw(iwm,ix,iwa)+dxstep 1106 if(abs(dxstep).gt.dxwmax) dxwmax=abs(dxstep) 1107 endif 1108 46 continue 1109 45 continue 1110 44 continue 1111 if(dxwmax.gt.dxmax) dxmax=dxwmax 1112 endif 1113 if(nsaloc.gt.0) then 1114 do 47 ix=1,3 1115 do 48 isa=1,nsaloc 1116 ys(isa,ix)=xs(isa,ix) 1117 dxstep=beta5*fs(isa,ix)/ws(isa) 1118 if(iand(isdt(isa),mfixed).ne.lfixed) then 1119 xs(isa,ix)=ys(isa,ix)+dxstep 1120 if(abs(dxstep).gt.dxsmax) dxsmax=abs(dxstep) 1121 endif 1122 48 continue 1123 47 continue 1124 if(dxsmax.gt.dxmax) dxmax=dxsmax 1125 endif 1126 call ga_dgop(mrg_d12,dxmax,1,'max') 1127c 1128c shake 1129c 1130 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 1131 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 1132c 1133 fnorm2=zero 1134c 1135 if(nwmloc.gt.0) then 1136 do 52 iwa=1,nwa 1137 do 53 ix=1,3 1138 do 54 iwm=1,nwmloc 1139 fw(iwm,ix,iwa)=(xw(iwm,ix,iwa)-yw(iwm,ix,iwa))*ww(iwa)/beta5 1140 pcgw(iwm,ix,iwa)=(yw(iwm,ix,iwa)-vw(iwm,ix,iwa))*ww(iwa)/beta4 1141 vw(iwm,ix,iwa)=yw(iwm,ix,iwa) 1142 fnorm2=fnorm2+fw(iwm,ix,iwa)**2 1143 54 continue 1144 53 continue 1145 52 continue 1146 endif 1147 if(nsaloc.gt.0) then 1148 do 55 ix=1,3 1149 do 56 isa=1,nsaloc 1150 fs(isa,ix)=(xs(isa,ix)-ys(isa,ix))*ws(isa)/beta5 1151 pcgs(isa,ix)=(ys(isa,ix)-vs(isa,ix))*ws(isa)/beta4 1152 vs(isa,ix)=ys(isa,ix) 1153 fnorm2=fnorm2+fs(isa,ix)**2 1154 56 continue 1155 55 continue 1156 endif 1157c 1158c global sum fnorm2 1159c 1160 call ga_dgop(mrg_d13,fnorm2,1,'+') 1161c 1162 beta1=sqrt(fnorm2/fnorm1) 1163 fnorm1=fnorm2 1164c 1165 ecg2=epot 1166 epcgw=epotw 1167 epcgsw=epotsw 1168 epcgs=epots 1169c 1170 ecgdif=ecg2-ecgold 1171 ecgold=ecg2 1172 edt=ecg2-ecg0 1173 ecg1=ecg2 1174 ecgnew=ecg2 1175c 1176c record to lfntrj 1177c 1178 lxw=frequency(mdstep,nfcoor) 1179 lxs=frequency(mdstep,nfscoo) 1180c 1181 if(lxw.or.lxs) then 1182 if(lqmd) then 1183 call qmd_wrttrj(lfntrj,mwm,nwmloc,mwa,nwa,msa,nsaloc, 1184 + .true.,.false.,.true.,.false.,box,stime,pres,temp,tempw,temps, 1185 + dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xs),dbl_mb(i_vs)) 1186 else 1187 call sp_wrttrj(lfntrj,lxw,.false.,.false.,lxs,.false.,.false., 1188 + stime,pres,temp,tempw,temps, 1189 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw), 1190 + dbl_mb(i_xwcr),int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs), 1191 + dbl_mb(i_fs)) 1192 endif 1193 endif 1194c 1195c write restart file 1196c 1197 ldone=.not.(icgit.lt.mcgit.and.dxmax.gt.dxcgmx.and.ecgdif.lt.zero) 1198c 1199 cqrs=' ' 1200 if(frequency(icgit,nfqrs).or.(ldone.and.eqrs.gt.epot)) then 1201 write(projct,4000) nserie,msdit,icgit,filnam(1:56) 1202 4000 format(i2,' em ',i7,' + ',i7,' ',a) 1203 call md_wrtrst(lfnqrs,filqrs,.false.) 1204 cqrs='*' 1205 eqrs=ecg1 1206 endif 1207 if(me.eq.0.and.frequency(icgit,nfprop)) call prp_record() 1208c 1209c print minimization step data to output 1210c 1211 if(me.eq.0) then 1212 write(lfnout,600) icgit,cqrs,ecgdif,ecg1,epcgw,epcgsw,epcgs,dxmax 1213 600 format(i7,1x,a1,3x,5(1pe13.5),0pf12.8) 1214 endif 1215c 1216 if(itest.gt.0) call md_test() 1217c 1218 call timer_stop(201) 1219 if(.not.ldone) goto 100 1220c 1221 ecgdif=edt 1222c 1223 return 1224 end 1225 subroutine md_md() 1226c 1227 implicit none 1228c 1229#include "md_common.fh" 1230#include "mafdecls.fh" 1231#include "global.fh" 1232c 1233 logical frequency 1234 real*8 timer_wall_total 1235 external frequency,timer_wall_total 1236c 1237 character*1 mdt 1238 logical ltmp 1239c 1240 lfirst=.true. 1241 lxw=.false. 1242 lvw=.false. 1243 lxs=.false. 1244 lvs=.false. 1245 lesp=.false. 1246c 1247c equilibration 1248c 1249 ltmp=lhop 1250 lhop=.false. 1251 lequi=.true. 1252 lprpmf=.false. 1253 do 1 iequi=kequi+1,mequi 1254c 1255 mdstep=mdstep+1 1256 stime=stime+tstep 1257 lpmfc=npmf.gt.1.and.iequi.gt.npmf 1258 call timer_start(201) 1259 call md_newton() 1260 call prp_proper(mdstep,stime,eww,dbl_mb(i_esw), 1261 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 1262 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 1263 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsmp)) 1264 call md_server 1265 call timer_stop(201) 1266 1 continue 1267 lpmfc=.true. 1268c 1269c data gathering 1270c 1271 lhop=ltmp 1272 lequi=.false. 1273 lprpmf=iprpmf.ne.0 1274 if(lprpmf) lfnpmf=-iabs(lfnpmf) 1275 call timer_start(205) 1276 do 2 idacq=kdacq+1,mdacq 1277c 1278 mdstep=mdstep+1 1279 stime=stime+tstep 1280c 1281 lxw=frequency(mdstep,nfcoor) 1282 lvw=frequency(mdstep,nfvelo) 1283 lfw=frequency(mdstep,nfforc) 1284 lxs=frequency(mdstep,nfscoo) 1285 lvs=frequency(mdstep,nfsvel) 1286 lfs=frequency(mdstep,nfsfor) 1287 lesp=frequency(mdstep,nfesp) 1288c 1289 call md_timer_init() 1290c 1291 call timer_start(201) 1292c 1293 call md_newton() 1294c 1295 call timer_start(6) 1296c 1297 if(lfw.or.lfs) then 1298 call sp_gaputf(me,dbl_mb(i_fw),nwmloc,dbl_mb(i_fs),nsaloc) 1299 endif 1300c 1301 call timer_stop(6) 1302c 1303 call timer_start(55) 1304 mdt=' ' 1305 if(iguide.gt.0) mdt='g' 1306 write(projct,4000) nserie,mdt,mequi,idacq,tmpext,prsext, 1307 + filnam(1:38) 1308 4000 format(i2,' md',a1,i7,' + ',i7,' @',f7.2,e9.2,' ',a) 1309 if(frequency(mdstep,nfrest)) then 1310 call md_wrtrst(lfnrst,rfile,.true.) 1311 endif 1312 if(frequency(mdstep,nftime)) call md_wrtime 1313c 1314 call timer_stop(55) 1315c 1316 call prp_proper(mdstep,stime,eww,dbl_mb(i_esw), 1317 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 1318 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 1319 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsmp)) 1320 call md_server 1321 call prp_step(mdstep,stime,eww,dbl_mb(i_esw), 1322 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm) 1323c 1324 call timer_stop(201) 1325c 1326 if(frequency(mdstep,nfnewf).and.idacq.ne.mdacq) then 1327 if(me.eq.0) call md_fopen(.true.) 1328 endif 1329c 1330c 1331 if(lstop.and.(.not.lqmmm)) then 1332 if(me.eq.0) write(*,1000) tleft,tneed 1333 1000 format(///,' Time left (',f12.3,' s) is less than twice the ', 1334 + ' time needed to reach writing the next restart file (', 1335 + f12.3,' s)',//,' Simulation aborted',//) 1336 return 1337 endif 1338c 1339 2 continue 1340c 1341 return 1342 end 1343 subroutine md_ti() 1344c 1345 implicit none 1346c 1347#include "md_common.fh" 1348#include "mafdecls.fh" 1349#include "global.fh" 1350#include "msgids.fh" 1351c 1352 logical frequency,prp_mcti_step,prp_rdmri,sp_rdmri 1353 integer prp_dfr 1354 external frequency,prp_mcti_step,prp_rdmri,sp_rdmri,prp_dfr 1355c 1356 logical done 1357 integer npp,i,irunp,jequi,jdacq,ndec,ndum 1358 real*8 rdum 1359 character*256 filrun 1360c 1361 if(nserie.eq.0) then 1362 if(mropt.ge.2) then 1363 if(me.eq.0) then 1364 open(unit=lfnmri,file=filmri(1:index(filmri,' ')-1), 1365 + form='unformatted',status='unknown',err=9999) 1366 read(lfnmri) npp 1367 if(npp.ne.np) call md_abort('Number of nodes changed',npp) 1368 endif 1369 krun=0 1370 endif 1371 if(me.eq.0) then 1372 open(unit=lfnmro,file=filmro(1:index(filmro,' ')-1), 1373 + form='unformatted',status='unknown',err=9999) 1374 write(lfnmro) np,mrun,mequi,mdacq 1375 endif 1376 else 1377c 1378 if(me.eq.0) then 1379 open(unit=lfnmro,file=filmro(1:index(filmro,' ')-1), 1380 + form='unformatted',status='unknown',err=9999) 1381 read(lfnmro) npp,mrun,mequi,mdacq 1382 if(npp.ne.np) call md_abort('Number of nodes changed',npp) 1383 endif 1384 krun=0 1385 do 11 irun=1,mrun 1386 if(me.eq.0) then 1387 read(lfnmro,end=12,err=12) irunp,kequi,kdacq 1388 endif 1389 if(.not.prp_rdmri(lfnmro,ndec,mropt)) goto 12 1390 if(.not.sp_rdmri(lfnmro,stime,pres,temp,tempw,temps, 1391 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 1392 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs))) goto 12 1393 krun=irun 1394 11 continue 1395 12 continue 1396c 1397 if(me.eq.0) then 1398 open(unit=lfnmri,file=filmri(1:index(filmri,' ')-1), 1399 + form='unformatted',status='unknown',err=9999) 1400 read(lfnmri) npp 1401 if(npp.ne.np) call md_abort('Number of nodes changed',npp) 1402 rewind(lfnmro) 1403 read(lfnmro) npp,mrun,mequi,mdacq 1404 if(npp.ne.np) call md_abort('Number of nodes changed',npp) 1405 rewind(lfngib) 1406 endif 1407c 1408 do 13 irun=1,mrun 1409 if(me.eq.0) then 1410 read(lfnmri) irunp,kequi,kdacq 1411 endif 1412 if(.not.prp_rdmri(lfnmri,ndec,mropt)) 1413 + call md_abort('Error in mri file',0) 1414 if(.not.sp_rdmri(lfnmri,stime,pres,temp,tempw,temps, 1415 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 1416 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs))) 1417 + call md_abort('Error in mri file',0) 1418 read(lfnmro) irunp,kequi,kdacq 1419 if(.not.prp_rdmri(lfnmro,ndec,mropt)) 1420 + call md_abort('Error in mro file',0) 1421 if(.not.sp_rdmri(lfnmro,stime,pres,temp,tempw,temps, 1422 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 1423 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs))) 1424 + call md_abort('Error in mro file',0) 1425 if(me.eq.0) then 1426 read(lfngib,2000) ndum 1427 2000 format(7x,i7) 1428 read(lfngib,2001) (rdum,i=1,24) 1429 2001 format(4e20.12) 1430 read(lfngib,2001) (rdum,i=1,ndum) 1431 read(lfngib,2001) (rdum,i=1,ndum) 1432 read(lfngib,2002) ndum 1433 2002 format(i10) 1434 read(lfngib,2003) rdum 1435 read(lfngib,2003) rdum 1436 2003 format(e20.12) 1437 if(ndec.gt.0) read(lfngib,2001) (rdum,i=1,nsa) 1438 endif 1439 13 continue 1440c 1441 endif 1442c 1443 call ga_brdcst(mrg_d30,krun,ma_sizeof(mt_int,1,mt_byte),0) 1444c 1445 do 1 irun=krun+1,mrun 1446c 1447 if(npg.gt.1) then 1448 if(irun.ne.ipg+1) goto 1 1449 endif 1450c 1451 if(me.eq.0) then 1452 if(nfcoor.gt.0.or.nfscoo.gt.0.or.nfvelo.gt.0.or.nfsvel.gt.0) then 1453 write(filrun,'(a,a,i5.5,a)') filtrj(1:index(filtrj,'.trj')-1),'-', 1454 + irun,'.trj ' 1455 open(unit=lfntrj,file=filrun(1:index(filrun,' ')-1), 1456 + form='formatted',status='unknown') 1457 call cf_trjhdr(lfntrj) 1458 endif 1459 if(nfprop.gt.0) then 1460 write(filrun,'(a,a,i5.5,a)') filprp(1:index(filprp,'.prp')-1),'-', 1461 + irun,'.prp ' 1462 open(unit=lfnprp,file=filrun(1:index(filrun,' ')-1), 1463 + form='formatted',status='unknown') 1464 endif 1465 if(nfrdf.gt.0) then 1466 write(filrun,'(a,a,i5.5,a)') filrdf(1:index(filprp,'.prp')-1),'-', 1467 + irun,'.rdf ' 1468 open(unit=lfnrdf,file=filrun(1:index(filrun,' ')-1), 1469 + form='formatted',status='unknown') 1470 endif 1471 endif 1472c 1473 if(irun.eq.1.and.iand(ivopt,2).eq.2) nfgaus=ivreas 1474 if(irun.eq.maxlam.and.iand(ivopt,4).eq.4) nfgaus=ivreas 1475c 1476 lfirst=.true. 1477c 1478c initialize parameters 1479c 1480 call cf_lambda(lamtyp,irun,maxlam,elam,lfnout,lfnpmf, 1481 + rlambd,dlambd,filnam) 1482c 1483c property initialization 1484c 1485 call prp_init() 1486c 1487 if(mropt.ge.2) then 1488 if(me.eq.0) then 1489 read(lfnmri,end=399,err=399) irunp,kequi,kdacq 1490 if(irunp.ne.irun) call md_abort('Number of run changed',irunp) 1491 if(.not.prp_rdmri(lfnmri,ndec,mropt)) 1492 + call md_abort('Error in mri',0) 1493 if(kequi+kdacq.lt.mequi) then 1494 kequi=kequi+kdacq 1495 kdacq=0 1496 elseif(kequi.lt.mequi) then 1497 if(kdacq+kequi-mequi.gt.0) then 1498 kdacq=prp_dfr(kdacq+kequi-mequi) 1499 else 1500 kdacq=0 1501 endif 1502 kequi=mequi 1503 endif 1504 goto 398 1505 399 continue 1506 kequi=0 1507 kdacq=0 1508 398 continue 1509 endif 1510 call ga_brdcst(mrg_d31,kequi,ma_sizeof(mt_int,1,mt_byte),0) 1511 call ga_brdcst(mrg_d32,kdacq,ma_sizeof(mt_int,1,mt_byte),0) 1512c 1513c kequi will be 0 if no records could be read from lfnmri, and 1514c the coordinates and velocities will be used from the previous 1515c lambda run. 1516c if this is the first run there must be something wrong 1517c 1518 if(kequi.eq.0.and.kdacq.eq.0) then 1519 if(irun.eq.1) call md_abort('No records found on mri',me) 1520 else 1521 if(.not.sp_rdmri(lfnmri,stime,pres,temp,tempw,temps, 1522 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 1523 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs))) 1524 + call md_abort('Error reading mri in sp_rdmri',0) 1525 endif 1526 endif 1527c 1528 lxw=.false. 1529 lvw=.false. 1530 lxs=.false. 1531 lvs=.false. 1532 lesp=.false. 1533c 1534 if(mropt.eq.2) then 1535 kequi=0 1536 kdacq=0 1537 call prp_init 1538 endif 1539c 1540c equilibration 1541c 1542 lequi=.true. 1543 jequi=kequi 1544 lprpmf=.false. 1545 iprpmf=0 1546 do 2 iequi=kequi+1,mequi 1547c 1548 mdstep=mdstep+1 1549 lpmfc=npmf.gt.1.and.iequi.gt.npmf 1550 call timer_start(201) 1551 call md_newton() 1552 call timer_stop(201) 1553 jequi=iequi 1554 stime=stime+tstep 1555 2 continue 1556 lpmfc=.true. 1557c 1558c data gathering 1559c 1560 mdstep=kdacq 1561 if(kdacq.eq.0) stime=zero 1562c 1563 ndec=0 1564 if(npgdec.gt.1) call cf_dera_init() 1565c 1566 lequi=.false. 1567 jdacq=kdacq 1568 lprpmf=.true. 1569 iprpmf=-1 1570 do 3 idacq=kdacq+1,mdacq 1571c 1572 mdstep=mdstep+1 1573c 1574 lxw=frequency(mdstep,nfcoor) 1575 lvw=frequency(mdstep,nfvelo) 1576 lfw=frequency(mdstep,nfforc) 1577 lxs=frequency(mdstep,nfscoo) 1578 lvs=frequency(mdstep,nfsvel) 1579 lfs=frequency(mdstep,nfsfor) 1580c 1581 call md_timer_init() 1582c 1583 call timer_start(201) 1584c 1585 call md_newton() 1586c 1587 if(npgdec.gt.1) ndec=ndec+1 1588c 1589 call cf_mcti_kin(int_mb(i_is+(lsatt-1)*msa), 1590 + int_mb(i_is+(lsgan-1)*msa),dbl_mb(i_vs),nsaloc) 1591c 1592 done=prp_mcti_step(idacq,ldacq) 1593c 1594 call prp_proper(mdstep,stime,eww,dbl_mb(i_esw), 1595 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 1596 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 1597 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsmp)) 1598 call prp_step(mdstep,stime,eww,dbl_mb(i_esw), 1599 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm) 1600c 1601 if(lfw.or.lfs) then 1602 call sp_gaputf(me,dbl_mb(i_fw),nwmloc,dbl_mb(i_fs),nsaloc) 1603 endif 1604c 1605 write(projct,4000) nserie,irun,mrun,mequi,idacq,tmpext, 1606 + filnam(1:32) 1607 4000 format(i2,' ti ',i5,'/',i5,' :',i7,'+ ',i7,' @ ',f7.2,' K ',a) 1608 if(frequency(mdstep,nfrest)) call md_wrtrst(lfnrst,rfile,.true.) 1609 if(frequency(mdstep,nftime)) call md_wrtime 1610c 1611 call timer_stop(201) 1612c 1613 jdacq=idacq 1614 if(idacq.gt.ldacq.and.done) goto 4 1615c 1616 stime=stime+tstep 1617c 1618 3 continue 1619 4 continue 1620c 1621 if(me.eq.0) then 1622 write(lfnmro) irun,jequi,jdacq 1623 endif 1624 call prp_wrtmro(lfnmro,ndec) 1625 if(me.eq.0) then 1626 call sp_wrtmro(lfnmro,stime,pres,temp,tempw,temps, 1627 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 1628 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),projct) 1629 endif 1630c 1631 call prp_mcti_run(rlambd,dlambd,ndec) 1632c 1633 if(ndec.gt.2) call cf_print_dera(lfnout,ndec) 1634c 1635 if(me.eq.0) then 1636 if(nfcoor.gt.0.or.nfscoo.gt.0.or.nfvelo.gt.0.or.nfsvel.gt.0) then 1637 close(unit=lfntrj,status='keep') 1638 endif 1639 if(nfprop.gt.0) then 1640 close(unit=lfnprp,status='keep') 1641 endif 1642 endif 1643c 1644 1 continue 1645c 1646 call prp_mcti(npgdec,filnam) 1647c 1648 if(me.eq.0) then 1649 if(nserie.eq.0) then 1650 if(mropt.ge.2) then 1651 close(unit=lfnmri,status='keep') 1652 endif 1653 close(unit=lfnmro,status='keep') 1654 else 1655 close(unit=lfnmro,status='keep') 1656 endif 1657 endif 1658c 1659 return 1660c 1661 9999 continue 1662 call md_abort('Failed to open file mro',0) 1663 return 1664 end 1665 subroutine md_newton() 1666c 1667 implicit none 1668c 1669#include "md_common.fh" 1670#include "mafdecls.fh" 1671c 1672 logical frequency 1673 external frequency 1674c 1675 external timer_wall,timer_wall_average,timer_wall_minimum 1676 real*8 timer_wall,timer_wall_average,timer_wall_minimum 1677 external timer_wall_total 1678 real*8 timer_wall_total 1679c 1680 real*8 dxmax,wallt,syntim 1681c 1682c start timer 1683c 1684 call timer_start(1) 1685c 1686c debug code 1687c 1688 if(idebug.gt.0) then 1689 write(lfndbg,'(a,i5,f12.6)') 'Timestep ',mdstep,stime 1690 endif 1691c 1692c end debug code 1693c 1694 if(tann2.gt.0.and.tmpext1.ne.tmpext2.and.tann1.lt.tann2)then 1695 if(stime.gt.tann1.and.stime.lt.tann2) then 1696 tmpext=((tann2-stime)*tmpext1+(stime-tann1)*tmpext2)/(tann2-tann1) 1697 else 1698 tmpext=tmpext1 1699 if(stime.ge.tann2) tmpext=tmpext2 1700 endif 1701 endif 1702c 1703 wallt=timer_wall(203) 1704c 1705c dynamic load balancing data 1706c 1707 if(itload.eq.0) then 1708 wallt=timer_wall(203) 1709 syntim=timer_wall(204) 1710 elseif(itload.eq.1) then 1711 wallt=timer_wall_minimum(203) 1712 syntim=timer_wall_minimum(204) 1713 elseif(itload.eq.2) then 1714 wallt=timer_wall_average(203) 1715 syntim=timer_wall_average(204) 1716 else 1717 wallt=timer_wall_average(203) 1718 syntim=timer_wall_average(204) 1719 if(me.eq.0) then 1720 wallt=timer_wall_minimum(203) 1721 syntim=timer_wall_minimum(204) 1722 endif 1723 endif 1724c 1725 if(me.eq.0.and.ioload.eq.1) then 1726 if(corrio.ge.syntim) then 1727 syntim=zero 1728 else 1729 syntim=syntim-corrio 1730 endif 1731 endif 1732c 1733 if(lpair) call timer_reset(203) 1734 call timer_start(203) 1735c 1736c reassign velocities 1737c 1738 if(frequency(mdstep,nfgaus)) then 1739 call cf_gauss(tgauss,frgaus,nwmloc,nsaloc, 1740 + dbl_mb(i_vw),dbl_mb(i_vs),int_mb(i_iw+(lwdyn-1)*mwm), 1741 + int_mb(i_is+(lsdyn-1)*msa),int_mb(i_is+(lsatt-1)*msa)) 1742 endif 1743c 1744 lpair=lfirst.or.frequency(mdstep,nfpair) 1745 lload=lfirst.or.frequency(mdstep,nfload) 1746 lhop=frequency(mdstep,nfhop) 1747 llong=(lfirst.or.frequency(mdstep,nflong).or.lpair).and.ltwin 1748c 1749 call timer_stop(1) 1750c 1751 if(lpair) then 1752c 1753 call timer_start(2) 1754c 1755c center of mass coordinates 1756c 1757 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 1758 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 1759 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 1760c 1761c periodic boundary conditions 1762c 1763 call md_fold(int_mb(i_iw),int_mb(i_is), 1764 + dbl_mb(i_xw),dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_xsm)) 1765c 1766 call timer_stop(2) 1767c 1768 if(lload) then 1769 if(.not.lqmd) then 1770c 1771 call timer_start(3) 1772c 1773 if(me.eq.0.and.ioload.eq.2) then 1774 if(corrio.ge.syntim.and. 1775 + ((nfcoor.gt.nfpair.and.nfcoor.gt.np).or. 1776 + (nfcoor.eq.0.and.nfscoo.gt.nfpair.and.nfscoo.gt.np))) then 1777 syntim=zero 1778 else 1779 syntim=syntim-corrio 1780 endif 1781 endif 1782c 1783 call sp_balanc(stime,syntim,wallt,frequency(mdstep,nfsync)) 1784 call timer_reset(204) 1785c 1786 call timer_stop(3) 1787 call timer_start(4) 1788c 1789c atom redistribution 1790c 1791 call sp_travel(box,dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 1792 + dbl_mb(i_gw),int_mb(i_iw),nwmloc,dbl_mb(i_xs),dbl_mb(i_vs), 1793 + dbl_mb(i_gs),int_mb(i_is),nsaloc) 1794c 1795 call timer_stop(4) 1796c 1797 endif 1798 endif 1799c 1800 endif 1801c 1802 call timer_start(5) 1803c 1804c center of mass coordinates 1805c 1806 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 1807 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 1808 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 1809c 1810c subtract center of mass coordinates from reference coordinates 1811c 1812 if(idifco.gt.0) 1813 + call md_addref(.false.,dbl_mb(i_xwm),dbl_mb(i_xwcr), 1814 + dbl_mb(i_xsm),dbl_mb(i_xscr),dbl_mb(i_dsr)) 1815c 1816 call timer_stop(5) 1817 call timer_start(6) 1818c 1819 if(lfw.or.lfs) then 1820 call sp_gaputf(me,dbl_mb(i_fw),nwmloc,dbl_mb(i_fs),nsaloc) 1821 call sp_wrttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs, 1822 + stime,pres,temp,tempw,temps, 1823 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw), 1824 + dbl_mb(i_xwcr),int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs), 1825 + dbl_mb(i_fs)) 1826 endif 1827c 1828 call timer_stop(6) 1829c 1830 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 1831 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 1832 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 1833c 1834 corrio=costio 1835 if(me.eq.0.and..not.lequi.and.(lxw.or.lvw.or.lxs.or.lvs)) then 1836c 1837 call timer_start(6) 1838c 1839 if(lqmd) then 1840 call qmd_wrttrj(lfntrj,mwm,nwmloc,mwa,nwa,msa,nsaloc, 1841 + .true.,.false.,.true.,.false.,box,stime-tstep,pres,temp,tempw, 1842 + temps, 1843 + dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xs),dbl_mb(i_vs)) 1844 else 1845 if(.not.lfw.and..not.lfs) then 1846 call sp_wrttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs, 1847 + stime,pres,temp,tempw,temps, 1848 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_fw), 1849 + dbl_mb(i_xwcr),int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs), 1850 + dbl_mb(i_fs)) 1851 endif 1852 endif 1853 call timer_stop(6) 1854c 1855 costio=max(costio,timer_wall(6)) 1856 corrio=costio-timer_wall(6) 1857c 1858 endif 1859c 1860c atomic forces and potential energies 1861c 1862 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 1863 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 1864c 1865c self-guided forces 1866c 1867 if(iguide.gt.0) then 1868 call timer_start(48) 1869 call md_guided(dbl_mb(i_fw),dbl_mb(i_fs), 1870 + dbl_mb(i_gw),dbl_mb(i_gs)) 1871 call timer_stop(48) 1872 endif 1873c 1874c center of mass options 1875c 1876 if(icmopt.gt.0) then 1877 call timer_start(48) 1878 call md_cmopt(dbl_mb(i_vs),dbl_mb(i_fs),dbl_mb(i_fcm), 1879 + int_mb(i_is+(lsmol-1)*msa),int_mb(i_is+(lsatt-1)*msa)) 1880 call timer_stop(48) 1881 endif 1882c 1883 if(imembr.ne.0) call md_membrane_forces(int_mb(i_mm), 1884 + dbl_mb(i_fm),dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_fs), 1885 + dbl_mb(i_wws)) 1886c 1887c time step 1888c 1889 call timer_start(49) 1890 call cf_mdstep(int_mb(i_iw+(lwdyn-1)*mwm),dbl_mb(i_xw), 1891 + dbl_mb(i_yw),dbl_mb(i_vw),dbl_mb(i_vwt),dbl_mb(i_fw),nwmloc, 1892 + int_mb(i_is+(lsdyn-1)*msa),int_mb(i_is+(lsatt-1)*msa), 1893 + dbl_mb(i_xs),dbl_mb(i_ys),dbl_mb(i_vs),dbl_mb(i_vst), 1894 + dbl_mb(i_fs),nsaloc,int_mb(i_iw+(lwgmn-1)*mwm), 1895 + int_mb(i_is+(lsgan-1)*msa),int_mb(i_is+(lssgm-1)*msa),tmpext, 1896 + int_mb(i_is+(lshop-1)*msa)) 1897 call timer_stop(49) 1898c 1899c shake 1900c 1901 call timer_start(50) 1902 call md_shake(dbl_mb(i_xw),dbl_mb(i_yw),int_mb(i_iw), 1903 + dbl_mb(i_xs),dbl_mb(i_ys),int_mb(i_is),dxmax) 1904 call timer_stop(50) 1905c 1906c recalculate velocities 1907c 1908 call cf_veloc(nwmloc,dbl_mb(i_xw),dbl_mb(i_yw),dbl_mb(i_vw), 1909 + nsaloc,dbl_mb(i_xs),dbl_mb(i_ys),dbl_mb(i_vs)) 1910c 1911c velocity scaling to preset temperature 1912c 1913 if(frequency(mdstep,nfgaus)) then 1914 call cf_vscale(tgauss,nwmloc,nsaloc, 1915 + dbl_mb(i_vw),dbl_mb(i_vwt),dbl_mb(i_vs),dbl_mb(i_vst), 1916 + int_mb(i_iw+(lwdyn-1)*mwm), 1917 + int_mb(i_is+(lsdyn-1)*msa),int_mb(i_is+(lsatt-1)*msa)) 1918 if(iand(ivopt,1).eq.1) nfgaus=0 1919 endif 1920c 1921c coordinate scaling 1922c 1923 call timer_start(51) 1924 call cf_final(dbl_mb(i_xw),dbl_mb(i_xwm),dbl_mb(i_yw), 1925 + dbl_mb(i_vw),dbl_mb(i_vwt),nwmloc, 1926 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_ys),dbl_mb(i_vs), 1927 + dbl_mb(i_vst),int_mb(i_is+(lsatt-1)*msa), 1928 + int_mb(i_is+(lsmol-1)*msa),int_mb(i_is+(lsdyn-1)*msa), 1929 + int_mb(i_is+(lsfrc-1)*msa),int_mb(i_is+(lshop-1)*msa), 1930 + dbl_mb(i_zs), 1931 + dbl_mb(i_esk),nsaloc,box,vlat,pres,temp,tempw,temps) 1932 call timer_stop(51) 1933c 1934c center of mass coordinates 1935c 1936 call timer_start(52) 1937c 1938 if(idifco.gt.0.or. 1939 + frequency(mdstep,nfcntr).or.frequency(mdstep,nfslow)) then 1940c 1941 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 1942 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 1943 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 1944c 1945c add center of mass coordinates from reference coordinates 1946c 1947 if(idifco.gt.0) 1948 + call md_addref(.true.,dbl_mb(i_xwm),dbl_mb(i_xwcr), 1949 + dbl_mb(i_xsm),dbl_mb(i_xscr),dbl_mb(i_dsr)) 1950c 1951c center solute in box 1952c 1953 if(frequency(mdstep,nfcntr)) then 1954 call cf_center(dbl_mb(i_xw),nwmloc,int_mb(i_is+(lsfrc-1)*msa), 1955 + dbl_mb(i_xs),nsaloc,idscb,nscb,icentr) 1956 endif 1957c 1958c remove overall translational motion 1959c 1960 if(frequency(mdstep,nfslow)) then 1961 call cf_slow(dbl_mb(i_xw),dbl_mb(i_vw),nwmloc,dbl_mb(i_xs), 1962 + dbl_mb(i_vs),int_mb(i_is+(lsatt-1)*msa),nsaloc) 1963 endif 1964c 1965 endif 1966 call timer_stop(52) 1967 call timer_start(53) 1968c 1969c update decomposition module 1970c 1971 if(.not.lqmd) call sp_update(me,vlat, 1972 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_xwcr),dbl_mb(i_vw),nwmloc, 1973 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),nsaloc) 1974c 1975 if(itest.gt.0) call md_test() 1976c 1977 lfirst=.false. 1978c 1979 call timer_stop(53) 1980 call timer_stop(203) 1981c 1982 return 1983 end 1984 subroutine md_addref(ladd,xwm,xwcr,xsm,xscr,dsr) 1985c 1986 implicit none 1987c 1988#include "md_common.fh" 1989c 1990 logical ladd 1991 real*8 xwm(mwm,3),xwcr(mwm,3),xsm(msm,3),xscr(msm,3),dsr(msm) 1992c 1993 integer i 1994c 1995c return 1996 if(ladd) then 1997 dwr=zero 1998 do 1 i=1,nwmloc 1999 xwcr(i,1)=xwcr(i,1)+xwm(i,1) 2000 xwcr(i,2)=xwcr(i,2)+xwm(i,2) 2001 xwcr(i,3)=xwcr(i,3)+xwm(i,3) 2002 dwr=dwr+xwcr(i,1)*xwcr(i,1)+xwcr(i,2)*xwcr(i,2)+ 2003 + xwcr(i,3)*xwcr(i,3) 2004 1 continue 2005 do 2 i=1,nsm 2006 xscr(i,1)=xscr(i,1)+xsm(i,1) 2007 xscr(i,2)=xscr(i,2)+xsm(i,2) 2008 xscr(i,3)=xscr(i,3)+xsm(i,3) 2009 dsr(i)=xscr(i,1)*xscr(i,1)+xscr(i,2)*xscr(i,2)+xscr(i,3)*xscr(i,3) 2010 2 continue 2011 else 2012 do 3 i=1,nwmloc 2013 xwcr(i,1)=xwcr(i,1)-xwm(i,1) 2014 xwcr(i,2)=xwcr(i,2)-xwm(i,2) 2015 xwcr(i,3)=xwcr(i,3)-xwm(i,3) 2016 3 continue 2017 do 4 i=1,nsm 2018 xscr(i,1)=xscr(i,1)-xsm(i,1) 2019 xscr(i,2)=xscr(i,2)-xsm(i,2) 2020 xscr(i,3)=xscr(i,3)-xsm(i,3) 2021 4 continue 2022 endif 2023c 2024 return 2025 end 2026 subroutine md_shake(xw,yw,iwl,xs,ys,isl,dmax) 2027c 2028 implicit none 2029c 2030#include "md_common.fh" 2031#include "mafdecls.fh" 2032#include "msgids.fh" 2033#include "global.fh" 2034c 2035 logical cf_shakep 2036 external cf_shakep 2037c 2038 real*8 xw(mwm,3,mwa),yw(mwm,3,mwa),xs(msa,3),ys(msa,3) 2039 integer iwl(mwm,miw2),isl(msa,mis2) 2040 real*8 dmax 2041c 2042 integer i,j,lhandl,ibbl,iwfr,iwto,isfr,isto,nbbl 2043 logical lself 2044c 2045 dmax=zero 2046c 2047 if(nwmloc.gt.0) then 2048 call cf_shakew(xw,yw,iwl(1,lwgmn),iwl(1,lwdyn),nwmloc) 2049 do 1 j=1,nwa 2050 do 2 i=1,nwmloc 2051 dmax=max(dmax,(xw(i,1,j)-yw(i,1,j))**2+ 2052 + (xw(i,2,j)-yw(i,2,j))**2+(xw(i,3,j)-yw(i,3,j))**2) 2053 2 continue 2054 1 continue 2055 endif 2056c 2057 3 continue 2058 if(nsaloc.gt.0) then 2059 call sp_nbbl(nbbl) 2060 do 4 ibbl=1,nbbl 2061 call sp_gethdl(ibbl,lhandl,lself,iwfr,iwto,isfr,isto) 2062 if(lself) call cf_shakes(lhandl,xs,ys,isl(1,lsgan),isl(1,lsatt), 2063 + isl(1,lssgm),isl(1,lsdyn),isl(1,lshop),isfr,isto) 2064 4 continue 2065 do 5 i=1,nsaloc 2066 dmax=max(dmax,(xs(i,1)-ys(i,1))**2+(xs(i,2)-ys(i,2))**2+ 2067 + (xs(i,3)-ys(i,3))**2) 2068 5 continue 2069 endif 2070 if(npmf.eq.1.or.lpmfc) then 2071 if(.not.cf_shakep(xs,ys,isl(1,lsgan),isl(1,lsatt),isl(1,lsdyn), 2072 + isl(1,lshop),nsaloc)) goto 3 2073 endif 2074c 2075 dmax=sqrt(dmax) 2076 call ga_dgop(mrg_d45,dmax,1,'max') 2077c 2078 return 2079 end 2080 subroutine md_fold(iwl,isl,xw,xwm,xs,xsm) 2081c 2082 implicit none 2083c 2084#include "md_common.fh" 2085c 2086 integer iwl(mwm,miw2),isl(msa,mis2) 2087 real*8 xw(mwm,3,mwa),xwm(mwm,3),xs(msa,3),xsm(msm,3) 2088c 2089 call cf_fold(nwmloc,xw,xwm,nsaloc,isl(1,lsatt),isl(1,lsmol), 2090 + xs,xsm) 2091c 2092 return 2093 end 2094 subroutine md_finit(iwl,isl,xw,xwm,xs,fw,fs,xsm,xsmp) 2095c 2096 implicit none 2097c 2098#include "md_common.fh" 2099#include "mafdecls.fh" 2100#include "global.fh" 2101c 2102 integer iwl(mwm,miw2),isl(msa,mis2) 2103 real*8 xw(mwm,3,mwa),xwm(mwm,3),xs(msa,3) 2104 real*8 fw(mwm,3,mwa,2),fs(msa,3,2) 2105 real*8 xsm(msm,3),xsmp(msm,3) 2106c 2107 integer i 2108c 2109 call timer_start(7) 2110c 2111 do 1 i=1,nsm 2112 xsmp(i,1)=xsm(i,1) 2113 xsmp(i,2)=xsm(i,2) 2114 xsmp(i,3)=xsm(i,3) 2115 1 continue 2116c 2117c initialize cafe 2118c 2119 call cf_init(stime,lpair,llong,box,vlat,vlati,zw,dbl_mb(i_zs), 2120 + eww,dbl_mb(i_esw),dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esa)) 2121c 2122 if(lpair) call md_zinit(int_mb(i_iwz),int_mb(i_isz)) 2123c 2124 call timer_stop(7) 2125c 2126 if(lqmd) then 2127 call timer_start(8) 2128 call qmd_forces(irtdb,int_mb(i_is+(lsgan-1)*msa), 2129 + int_mb(i_is+(lsatt-1)*msa),dbl_mb(i_xs),dbl_mb(i_fs),msa, 2130 + nsaloc,uqmd,lesp) 2131 uqmd=uqmd-uqmatm 2132 call timer_stop(8) 2133 return 2134 endif 2135c 2136 call timer_start(9) 2137 call ga_sync() 2138 call timer_stop(9) 2139c 2140 call timer_start(10) 2141c 2142 call sp_initf(fw,fs,llong,int_mb(i_iwz),int_mb(i_isz),lpair) 2143c 2144 call sp_putix(me,iwl,xw,nwmloc,isl,xs,nsaloc) 2145c 2146 call timer_stop(10) 2147c 2148 call timer_start(11) 2149 call ga_sync() 2150 call timer_stop(11) 2151c 2152 if(ncoll.gt.0) then 2153 call timer_start(12) 2154 call cf_collapse(ncoll,fcoll,nsaloc,nwmloc,isl(1,lsmol), 2155 + isl(1,lssgm), 2156 + dbl_mb(i_xs),dbl_mb(i_xsm),mst,dbl_mb(i_tsm),dbl_mb(i_fs), 2157 + dbl_mb(i_xw),dbl_mb(i_xwm),dbl_mb(i_fw)) 2158 call timer_stop(12) 2159 endif 2160 if(ifield.gt.0) then 2161 call timer_start(13) 2162 call cf_extern(stime,nsaloc,dbl_mb(i_fs), 2163 + isl(1,lsct1),nwmloc,dbl_mb(i_fw)) 2164 call timer_stop(13) 2165 endif 2166 call cf_multi(nsaloc,dbl_mb(i_xs),dbl_mb(i_fs),isl(1,lsgan), 2167 + isl(1,lsatt),isl(1,lsfrc),isl(1,lsdyn),isl(1,lsct1), 2168 + dbl_mb(i_ess),dbl_mb(i_fss),lfnpmf,lprpmf,iprpmf, 2169 + npmf.eq.1.or.lpmfc) 2170c 2171 if(ipolt.gt.0) then 2172 call timer_start(23) 2173 call md_induce(iwl,isl,xw,xwm,xs,dbl_mb(i_pw), 2174 + dbl_mb(i_pwp),dbl_mb(i_ps),dbl_mb(i_psp)) 2175 call timer_stop(23) 2176 endif 2177c 2178 if(lpme.and.llong) then 2179 if(lpert2) then 2180 call pme_energy(2,dbl_mb(i_xw),nwmloc,dbl_mb(i_xs), 2181 + isl(1,lsct1),isl(1,lssgm),nsaloc,epme(2)) 2182 call timer_start(24) 2183 call pme_init() 2184 call timer_stop(24) 2185 endif 2186 if(lpert3) then 2187 call pme_energy(3,dbl_mb(i_xw),nwmloc,dbl_mb(i_xs), 2188 + isl(1,lsct1),isl(1,lssgm),nsaloc,epme(3)) 2189 call timer_start(24) 2190 call pme_init() 2191 call timer_stop(24) 2192 endif 2193 call pme_chgrid(iset,dbl_mb(i_xw),nwmloc,dbl_mb(i_xs), 2194 + isl(1,lsct1),isl(1,lssgm),nsaloc,epme(iset)) 2195 endif 2196c 2197 return 2198 end 2199 subroutine md_forces(iwl,isl,xw,xwm,xs,fw,fs) 2200c 2201 implicit none 2202c 2203#include "md_common.fh" 2204#include "mafdecls.fh" 2205#include "global.fh" 2206c 2207 integer iwl(mwm,miw2),isl(msa,mis2) 2208 real*8 xw(mwm,3,mwa),xwm(mwm,3),xs(msa,3) 2209 real*8 fw(mwm,3,mwa,2),fs(msa,3,2) 2210c 2211 integer i,j,k 2212c 2213 if(.not.lqmd) then 2214c 2215 call md_fclass(iwl,isl,xw,xwm,xs,fw,fs) 2216c 2217 call timer_start(44) 2218 call timer_start(202) 2219 call timer_start(204) 2220 call ga_sync() 2221 call timer_stop(204) 2222 call timer_stop(202) 2223 call timer_stop(44) 2224c 2225 call timer_start(45) 2226 call sp_final(fw,fs,lpair,int_mb(i_iwz),int_mb(i_isz)) 2227 call timer_stop(45) 2228 2229 if(lqmmm) then 2230 call timer_start(46) 2231 call qmmm_forces(irtdb,mwm,nwmloc,mwa,nwa,int_mb(i_iwz),xw,fw, 2232 + msa,nsaloc,int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsdyn-1)*msa), 2233 + int_mb(i_is+(lsct1-1)*msa),int_mb(i_isz),xs,fs,uqmmm) 2234c uqmmm=uqmmm-uqmatm 2235 call timer_stop(46) 2236 endif 2237 2238 if(ltwin) then 2239 do 1 j=1,3 2240 do 2 k=1,nwa 2241 do 3 i=1,nwmloc 2242 fw(i,j,k,1)=fw(i,j,k,1)+fw(i,j,k,2) 2243 3 continue 2244 2 continue 2245 do 4 i=1,nsaloc 2246 fs(i,j,1)=fs(i,j,1)+fs(i,j,2) 2247 4 continue 2248 1 continue 2249 endif 2250c 2251 endif 2252c 2253 2254 call timer_start(47) 2255 call cf_fnorm(iwl(1,lwdyn),fw,nwmloc, 2256 + isl(1,lsatt),isl(1,lsdyn),fs,nsaloc,fnorm,fmax) 2257 call timer_stop(47) 2258c 2259 return 2260 end 2261 subroutine md_guided(fw,fs,gw,gs) 2262c 2263 implicit none 2264c 2265#include "md_common.fh" 2266#include "mafdecls.fh" 2267#include "global.fh" 2268c 2269 real*8 fw(mwm,3,mwa,2),fs(msa,3,2) 2270 real*8 gw(mwm,3,mwa),gs(msa,3) 2271c 2272 integer i,j,k 2273c 2274 do 1 k=1,nwa 2275 do 2 j=1,3 2276 do 3 i=1,nwmloc 2277 gw(i,j,k)=factgg*gw(i,j,k)+factgf*(fw(i,j,k,1)+fguide*gw(i,j,k)) 2278 fw(i,j,k,1)=fw(i,j,k,1)+fguide*gw(i,j,k) 2279 3 continue 2280 2 continue 2281 1 continue 2282c 2283 do 4 j=1,3 2284 do 5 i=1,nsaloc 2285 gs(i,j)=factgg*gs(i,j)+factgf*(fs(i,j,1)+fguide*gs(i,j)) 2286 fs(i,j,1)=fs(i,j,1)+fguide*gs(i,j) 2287 5 continue 2288 4 continue 2289c 2290 call sp_copyg(dbl_mb(i_gw),dbl_mb(i_gs)) 2291 call ga_sync() 2292c 2293 return 2294 end 2295 subroutine md_cmopt(vs,fs,fcm,ismol,isatt) 2296c 2297 implicit none 2298c 2299#include "md_common.fh" 2300#include "mafdecls.fh" 2301#include "global.fh" 2302c 2303 real*8 vs(msa,3),fs(msa,3,2),fcm(msm,5) 2304 integer ismol(msa),isatt(msa) 2305c 2306 call cf_scmfor(icmopt,ismol,isatt,vs,fs,nsaloc,fcm) 2307c 2308 return 2309 end 2310 subroutine md_zinit(iwz,isz) 2311c 2312 implicit none 2313c 2314#include "md_common.fh" 2315c 2316 integer iwz(mwm),isz(msa) 2317c 2318 integer i 2319c 2320 do 1 i=1,mwm 2321 iwz(i)=0 2322 1 continue 2323c 2324 do 2 i=1,msa 2325 isz(i)=0 2326 2 continue 2327c 2328 return 2329 end 2330 subroutine md_fclass(iwl,isl,xw,xwm,xs,fw,fs) 2331c 2332 implicit none 2333c 2334#include "md_common.fh" 2335#include "mafdecls.fh" 2336#include "global.fh" 2337c 2338 logical sp_local,cf_hopping 2339 external sp_local,cf_hopping 2340c 2341 integer iwl(mwm,miw2),isl(msa,mis2) 2342 real*8 xw(mwm,3,mwa),xwm(mwm,3),xs(msa,3) 2343 real*8 fw(mwm,3,mwa,2),fs(msa,3,2) 2344 logical lself,local,lpbcs 2345c 2346 integer ibbl,lhandl,iwfr,iwto,jwfr,jwto,isfr,isto,jsfr,jsto 2347 integer nbbl 2348 integer itcom1,itforc,itcom2 2349c 2350 logical lforce,ldone,lnew 2351c 2352 local=.true. 2353 itcom1=33 2354 itforc=34 2355 itcom2=35 2356c 2357 if(nbget.eq.0) then 2358 call sp_nbbl(nbbl) 2359 else 2360 itforc=35 2361 call timer_start(33) 2362 if(ipolt.eq.0) then 2363 if(nbget.lt.0) then 2364 call sp_prefetch_all(nbbl,iwl,xw,isl,xs) 2365 else 2366 call sp_prefetch(nbbl,iwl,xw,isl,xs) 2367 endif 2368 else 2369 call sp_prefetch_p(nbbl,iwl,xw,dbl_mb(i_pw),dbl_mb(i_pwp), 2370 + isl,xs,dbl_mb(i_ps),dbl_mb(i_psp)) 2371 endif 2372 call timer_stop(33) 2373 endif 2374c 2375 lforce=.true. 2376 ldone=.true. 2377c 2378 if(lpair.and.lhop) then 2379 lforce=.false. 2380 ldone=.false. 2381 endif 2382c 2383 100 continue 2384c 2385 do 1 ibbl=1,nbbl 2386c 2387 if(local.and..not.sp_local(ibbl)) then 2388 itcom1=36 2389 itforc=37 2390 itcom2=38 2391 local=.false. 2392 if(nbget.ne.0) itforc=36 2393 endif 2394c 2395 if(nbget.eq.0) then 2396 call timer_start(itcom1) 2397 if(ipolt.eq.0) then 2398 call sp_getxbl(ibbl,lhandl, 2399 + iwl,xw,iwfr,iwto,jwfr,jwto,isl,xs,isfr,isto,jsfr,jsto, 2400 + lself,lpbcs) 2401 else 2402 call sp_getxpbl(ibbl,lhandl, 2403 + iwl,xw,dbl_mb(i_pw),dbl_mb(i_pwp),iwfr,iwto,jwfr,jwto, 2404 + isl,xs,dbl_mb(i_ps),dbl_mb(i_psp),isfr,isto,jsfr,jsto, 2405 + lself,lpbcs) 2406 endif 2407 call timer_stop(itcom1) 2408 else 2409 call timer_start(34) 2410 call sp_nbwait(ibbl,lnew,lhandl,lself,lpbcs, 2411 + iwfr,iwto,jwfr,jwto,isfr,isto,jsfr,jsto,iwl,isl) 2412 call timer_stop(34) 2413 if(nbget.gt.0.and.lnew) then 2414 call timer_start(33) 2415 call sp_prefetch_next(iwl,xw,isl,xs) 2416 call timer_stop(33) 2417 endif 2418 endif 2419c 2420 call timer_start(itforc) 2421 call cf_comw(xw,xwm,jwfr,jwto) 2422c 2423 if(ipolt.eq.0) then 2424 call forces(lself,lpbcs,xw,xwm,fw,zw,dbl_mb(i_rtos),iwl(1,lwdyn), 2425 + int_mb(i_iwz),iwfr,iwto,jwfr,jwto,xs,dbl_mb(i_xsm),fs, 2426 + dbl_mb(i_zs),isl(1,lsgan),isl(1,lsatt),isl(1,lsdyn),isl(1,lsgrp), 2427 + isl(1,lsfrc),isl(1,lsmol),isl(1,lssss),isl(1,lsct1),isl(1,lsct2), 2428 + isl(1,lsct3),isl(1,lssgm),isl(1,lshop),int_mb(i_isz), 2429 + isfr,isto,jsfr,jsto,lpbc,lhandl, 2430 + .true.,eww,dbl_mb(i_esw),dbl_mb(i_ess),dbl_mb(i_fss), 2431 + dbl_mb(i_esa),int_mb(i_lseq),lforce) 2432 else 2433 call forcep(lself,lpbcs,xw,xwm,fw,dbl_mb(i_pw),dbl_mb(i_pwp),zw, 2434 + dbl_mb(i_rtos),iwl(1,lwdyn),int_mb(i_iwz),iwfr,iwto,jwfr,jwto,xs, 2435 + dbl_mb(i_xsm),fs,dbl_mb(i_ps),dbl_mb(i_psp),dbl_mb(i_zs), 2436 + isl(1,lsgan),isl(1,lsatt),isl(1,lsdyn),isl(1,lsgrp),isl(1,lsfrc), 2437 + isl(1,lsmol),isl(1,lssss),isl(1,lsct1),isl(1,lsct2),isl(1,lsct3), 2438 + isl(1,lssgm),isl(1,lshop),int_mb(i_isz), 2439 + isfr,isto,jsfr,jsto,lpbc,lhandl, 2440 + .true.,eww,dbl_mb(i_esw),dbl_mb(i_ess),dbl_mb(i_fss), 2441 + dbl_mb(i_esa),int_mb(i_lseq)) 2442 endif 2443 call timer_stop(itforc) 2444c 2445 if(nbget.eq.0) then 2446 call timer_start(itcom2) 2447 call sp_accfbl(ibbl,lhandl,fw,fs, 2448 + lpair,int_mb(i_iwz),int_mb(i_isz)) 2449 call timer_stop(itcom2) 2450 else 2451 call timer_start(37) 2452 call sp_nbaccfbl(ibbl,lhandl,fw,fs, 2453 + lpair,int_mb(i_iwz),int_mb(i_isz)) 2454 call timer_stop(37) 2455 endif 2456c 2457 1 continue 2458c 2459 if(nbget.ne.0) then 2460 call timer_start(38) 2461 call sp_nbwaitf() 2462 call timer_stop(38) 2463 endif 2464c 2465 if(.not.ldone) then 2466 if(cf_hopping(lpbc,lpbcs,stime,isl,isl(1,lssgm),isl(1,lsgan), 2467 + isl(1,lsct3),isl(1,lshop),xs,nsaloc,lfnhop)) then 2468 if(lpair) then 2469 lpair=.false. 2470 lload=.false. 2471 lforce=.true. 2472 ldone=.false. 2473 else 2474 lpair=.true. 2475 lload=.true. 2476 lforce=.true. 2477 ldone=.true. 2478 call sp_update_i(nsaloc,int_mb(i_is),nwmloc,int_mb(i_iw)) 2479 endif 2480 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 2481 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 2482 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 2483 goto 100 2484 endif 2485 endif 2486c 2487 if(lpme.and.llong) then 2488 call pme_forces(fw(1,1,1,2),nwmloc,fs(1,1,2), 2489 + isl(1,lsct1),isl(1,lssgm),nsaloc) 2490 endif 2491c 2492 return 2493 end 2494 subroutine md_induce(iwl,isl,xw,xwm,xs,pw,pwp,ps,psp) 2495c 2496 implicit none 2497c 2498#include "md_common.fh" 2499#include "mafdecls.fh" 2500#include "global.fh" 2501#include "msgids.fh" 2502c 2503 integer iwl(mwm,miw2),isl(msa,mis2) 2504 real*8 xw(mwm,3,mwa),xwm(mwm,3),xs(msa,3) 2505 real*8 pw(mwm,3,mwa,2),ps(msa,3,2) 2506 real*8 pwp(mwm,3,mwa,2,2),psp(msa,3,2,2) 2507c 2508 logical lself,lpbcs 2509 integer ibbl,nbbl,lhandl 2510 integer iwfr,iwto,jwfr,jwto,isfr,isto,jsfr,jsto 2511 integer i,j,k 2512 real*8 pmax 2513c 2514c initialize induced fields to zero for first order polarization 2515c 2516 if(ipolt.eq.1) then 2517 if(nwmloc.gt.0) then 2518 do 1 k=1,nwa 2519 do 2 j=1,3 2520 do 3 i=1,nwmloc 2521 pw(i,j,k,1)=zero 2522 pw(i,j,k,2)=zero 2523 3 continue 2524 2 continue 2525 1 continue 2526 endif 2527 if(nsaloc.gt.0) then 2528 do 4 j=1,3 2529 do 5 i=1,nsaloc 2530 ps(i,j,1)=zero 2531 ps(i,j,2)=zero 2532 5 continue 2533 4 continue 2534 endif 2535 if(lpert2.or.lpert3) then 2536 if(nwmloc.gt.0) then 2537 do 6 k=1,nwa 2538 do 7 j=1,3 2539 do 8 i=1,nwmloc 2540 pwp(i,j,k,1,1)=zero 2541 pwp(i,j,k,1,2)=zero 2542 pwp(i,j,k,2,1)=zero 2543 pwp(i,j,k,2,2)=zero 2544 8 continue 2545 7 continue 2546 6 continue 2547 endif 2548 if(nsaloc.gt.0) then 2549 do 9 j=1,3 2550 do 10 i=1,nsaloc 2551 psp(i,j,1,1)=zero 2552 psp(i,j,1,2)=zero 2553 psp(i,j,2,1)=zero 2554 psp(i,j,2,2)=zero 2555 10 continue 2556 9 continue 2557 endif 2558 endif 2559 endif 2560c 2561c iterative cycle to generate induced fields 2562c ------------------------------------------ 2563c 2564 npolit=0 2565 call sp_nbbl(nbbl) 2566 11 continue 2567 npolit=npolit+1 2568c 2569c copy fields from previous iteration 2570c ----------------------------------- 2571c 2572 if(nwmloc.gt.0) then 2573 do 12 k=1,nwa 2574 do 13 j=1,3 2575 do 14 i=1,nwmloc 2576 pw(i,j,k,2)=pw(i,j,k,1) 2577 pw(i,j,k,1)=zero 2578 14 continue 2579 13 continue 2580 12 continue 2581 endif 2582 if(nsaloc.gt.0) then 2583 do 15 j=1,3 2584 do 16 i=1,nsaloc 2585 ps(i,j,2)=ps(i,j,1) 2586 ps(i,j,1)=zero 2587 16 continue 2588 15 continue 2589 endif 2590 if(mdtype.gt.3) then 2591 if(nwmloc.gt.0) then 2592 do 17 k=1,nwa 2593 do 18 j=1,3 2594 do 19 i=1,nwmloc 2595 pwp(i,j,k,1,2)=pwp(i,j,k,1,1) 2596 pwp(i,j,k,1,1)=zero 2597 pwp(i,j,k,2,2)=pwp(i,j,k,2,1) 2598 pwp(i,j,k,2,1)=zero 2599 19 continue 2600 18 continue 2601 17 continue 2602 endif 2603 if(nsaloc.gt.0) then 2604 do 20 j=1,3 2605 do 21 i=1,nsaloc 2606 psp(i,j,1,2)=psp(i,j,1,1) 2607 psp(i,j,1,1)=zero 2608 psp(i,j,2,2)=psp(i,j,2,1) 2609 psp(i,j,2,1)=zero 2610 21 continue 2611 20 continue 2612 endif 2613 endif 2614c 2615c copy current fields into local global array 2616c ------------------------------------------- 2617c 2618 call sp_putp(me,pw,pwp,nwmloc,ps,psp,nsaloc,lpert2.or.lpert3) 2619c 2620c synchronize to ensure induced fields are available 2621c -------------------------------------------------- 2622c 2623 call ga_sync() 2624c 2625 do 22 ibbl=1,nbbl 2626c 2627 call sp_getxpbl(ibbl,lhandl, 2628 + iwl,xw,pw,pwp,iwfr,iwto,jwfr,jwto, 2629 + isl,xs,ps,psp,isfr,isto,jsfr,jsto,lself,lpbcs) 2630c 2631 call cf_comw(xw,xwm,jwfr,jwto) 2632c 2633 call induce(lself,lpbcs,xw,xwm,pw,pwp,iwl(1,lwdyn),int_mb(i_iwz), 2634 + iwfr,iwto,jwfr,jwto,xs,dbl_mb(i_xsm),ps,psp, 2635 + isl(1,lsgan),isl(1,lsatt),isl(1,lsdyn),isl(1,lsgrp),isl(1,lsfrc), 2636 + isl(1,lsmol),isl(1,lssss),isl(1,lsct1),isl(1,lsct2),isl(1,lsct3), 2637 + int_mb(i_isz),isfr,isto,jsfr,jsto,lpbc,lhandl) 2638c 2639 call sp_accpbl(ibbl,lhandl,pw,pwp,ps,psp, 2640 + lpair,int_mb(i_iwz),int_mb(i_isz)) 2641c 2642 22 continue 2643 lpair=.false. 2644 lload=.false. 2645c 2646 if(np.gt.0) call ga_sync() 2647c 2648 if(nwmloc.gt.0) then 2649 do 26 k=1,nwa 2650 do 27 j=1,3 2651 do 28 i=1,nwmloc 2652 pw(i,j,k,2)=pw(i,j,k,1) 2653 pw(i,j,k,1)=zero 2654 28 continue 2655 27 continue 2656 26 continue 2657 endif 2658 if(nsaloc.gt.0) then 2659 do 29 j=1,3 2660 do 30 i=1,nsaloc 2661 ps(i,j,2)=ps(i,j,1) 2662 ps(i,j,1)=zero 2663 30 continue 2664 29 continue 2665 endif 2666c 2667 call sp_getp(me,pw,pwp,nwmloc,ps,psp,nsaloc,lpert2.or.lpert3,1) 2668c 2669 if(nwmloc.gt.0) then 2670 do 31 k=1,nwa 2671 do 32 j=1,3 2672 do 33 i=1,nwmloc 2673 pw(i,j,k,1)=pw(i,j,k,1)+pw(i,j,k,2) 2674 pw(i,j,k,2)=zero 2675 33 continue 2676 32 continue 2677 31 continue 2678 endif 2679 if(nsaloc.gt.0) then 2680 do 34 j=1,3 2681 do 35 i=1,nsaloc 2682 ps(i,j,1)=ps(i,j,1)+ps(i,j,2) 2683 ps(i,j,2)=zero 2684 35 continue 2685 34 continue 2686 endif 2687c 2688 call sp_getp(me,pw,pwp,nwmloc,ps,psp,nsaloc,lpert2.or.lpert3,2) 2689c 2690 pmax=0.0d0 2691 do 23 k=1,nwa 2692 do 24 j=1,3 2693 do 25 i=1,nwmloc 2694 pmax=max(pmax,abs(pw(i,j,k,2)-pw(i,j,k,1))) 2695 25 continue 2696 24 continue 2697 23 continue 2698c 2699 if(np.gt.1) call ga_dgop(mrg_d06,pmax,1,'max') 2700c 2701 if(pmax.gt.ptol.and.npolit.le.mpolit.and.ipolt.gt.1) goto 11 2702c 2703c copy current fields into local global array 2704c ------------------------------------------- 2705c 2706 call sp_putp(me,pw,pwp,nwmloc,ps,psp,nsaloc,lpert2.or.lpert3) 2707c 2708 return 2709 end 2710 subroutine md_wrtrst(lfn,fil,lveloc) 2711c 2712 implicit none 2713c 2714#include "md_common.fh" 2715#include "mafdecls.fh" 2716#include "msgids.fh" 2717#include "util.fh" 2718c 2719 real*8 timer_wall 2720 external timer_wall 2721c 2722 integer lfn 2723 character*255 fil 2724 logical lveloc 2725 integer i,left 2726 character*255 filn 2727c 2728 if(me.eq.0) then 2729 if(keepr.eq.0) then 2730 open(unit=lfn,file=fil(1:index(fil,' ')-1), 2731 + form='formatted',status='unknown',err=9999) 2732 else 2733 if(ntype.eq.3.and.npg.le.1) then 2734 print*,npg,ntype,irun 2735 print*,fil(1:index(fil,' ')-1) 2736 print*,fil(1:index(fil,'.rst')-1) 2737 write(filn,'(a,i5.5,a)') fil(1:index(fil,'.rst')-1), 2738 + irun,'.rst ' 2739 else 2740 write(filn,'(a,a,i5.5,a)') fil(1:index(fil,'.rst')-1),'-', 2741 + keepr,'.rst ' 2742 endif 2743 open(unit=lfn,file=filn(1:index(filn,' ')-1), 2744 + form='formatted',status='unknown',err=9999) 2745 endif 2746 endif 2747c 2748 call sp_wrtrst(lfn,fil,lveloc,pres,temp,tempw,temps, 2749 + int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_gw), 2750 + dbl_mb(i_xwcr),int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs), 2751 + dbl_mb(i_gs),dbl_mb(i_xscr),projct,int_mb(i_lseq)) 2752c 2753 call md_wtrest(lfn) 2754c 2755 call prp_wtrest(lfn) 2756c 2757 call sp_wtrest(lfn) 2758c 2759 if(me.eq.0) then 2760 close(unit=lfn) 2761 endif 2762c 2763 if(keepr.gt.0) keepr=keepr+1 2764c 2765 call timer_stop(205) 2766 tneed=timer_wall(205) 2767 call timer_reset(205) 2768 call timer_start(205) 2769 left=util_time_remaining(irtdb) 2770 tleft=dble(left) 2771 i=0 2772 if(left.lt.0) then 2773 i=1 2774 elseif(tleft.gt.two*tneed) then 2775 i=1 2776 endif 2777 call ga_brdcst(mrg_d48,i,ma_sizeof(mt_int,1,mt_byte),0) 2778 lstop=i.eq.0 2779c 2780 return 2781c 2782 9999 continue 2783 call md_abort('Unable to open restart for writing',me) 2784 return 2785 end 2786 subroutine md_server 2787c 2788 implicit none 2789c 2790#include "md_common.fh" 2791c 2792#if !defined(WIN32) 2793 integer*4 create_client_socket 2794 integer client_socket_write 2795 external create_client_socket,client_socket_write 2796c 2797 integer*4 iiport,lens 2798c 2799 character*255 string 2800 integer numbyt 2801c 2802 if(iport.le.0.or.me.gt.0) return 2803c 2804c open socket to server 2805c 2806 iiport=iport 2807 if(.not.lserver) then 2808 write(*,1) server(1:index(server,' ')-1),iiport 2809 1 format('Attempt to open socket to ',a,' port ',i5) 2810 isocket=create_client_socket(server,iiport) 2811 lserver=isocket.gt.0 2812 endif 2813c 2814 if(lserver) then 2815 print*,'server socket open' 2816 else 2817 print*,'server socket error' 2818 endif 2819c 2820 if(.not.lserver) return 2821c 2822 write(string,1000) stime,etot,temp,pres*1.0125d-5,volume 2823 1000 format("tETpV",4f12.3,f12.6) 2824c 2825 lens=66 2826 string(66:66)=char(13) 2827c 2828 numbyt=client_socket_write(isocket,string,lens) 2829c 2830 print*,'Bytes written to socket is ',numbyt 2831c 2832#endif 2833 return 2834 end 2835 subroutine md_timer_init() 2836c 2837 implicit none 2838c 2839#include "md_common.fh" 2840c 2841 integer i 2842c 2843 do 1 i=1,200 2844 call timer_reset(i) 2845 1 continue 2846c 2847 return 2848 end 2849 subroutine md_wrtime 2850c 2851 implicit none 2852c 2853#include "md_common.fh" 2854#include "global.fh" 2855#include "msgids.fh" 2856c 2857 external timer_wall,timer_wall_total 2858 real*8 timer_wall,timer_wall_total 2859c 2860 integer i,j 2861 real*8 tim(56,1024) 2862c 2863 do 1 i=1,np 2864 do 2 j=1,56 2865 tim(j,i)=zero 2866 2 continue 2867 1 continue 2868c 2869 do 3 i=1,55 2870 tim(i,me+1)=timer_wall_total(i) 2871 tim(56,me+1)=tim(56,me+1)+tim(i,me+1) 2872 3 continue 2873c 2874 call ga_dgop(mrg_d04,tim,56*np,'+') 2875c 2876 if(me.eq.0) then 2877 write(lfntim,1000) stime 2878 1000 format('timings',/,f12.6) 2879 do 1002 j=1,np 2880 write(lfntim,1001) (tim(i,j),i=1,56) 2881 1001 format(10f7.3) 2882 1002 continue 2883 endif 2884c 2885 return 2886 end 2887 subroutine md_test 2888c 2889 implicit none 2890c 2891#include "md_common.fh" 2892c 2893 itest=itest-1 2894c 2895 if(me.ne.0) return 2896c 2897 if(iquant.eq.0) then 2898c 2899 if(ntype.eq.1) then 2900 write(lfntst,1000) etot 2901 1000 format('Energy = ',1pe12.3) 2902 endif 2903c 2904 if(ntype.eq.2) then 2905 write(lfntst,1100) stime,temp,volume,pres,etot 2906 1100 format('Time = ',0pf9.3,/, 2907 + 'Temperature = ',0pf8.2,/, 2908 + 'Volume = ',0pf8.2,/, 2909 + 'Pressure = ',1pe12.2,/, 2910 + 'Energy = ',1pe12.3,/) 2911 endif 2912c 2913 if(ntype.eq.3) write(lfntst,1200) isdit+icgit,etot 2914 1200 format('Iteration = ',i10,/, 2915 + 'Energy = ',1pe12.3,/) 2916c 2917 else 2918c 2919 if(ntype.eq.1) then 2920 write(lfntst,2000) etot 2921 2000 format('Energy = ',1pe12.3) 2922 endif 2923c 2924 if(ntype.eq.2) then 2925 write(lfntst,2100) stime,temp,etot 2926 2100 format('Time = ',0pf9.3,/, 2927 + 'Temperature = ',0pf8.2,/, 2928 + 'Energy = ',1pe12.3,/) 2929 endif 2930c 2931 endif 2932c 2933 if(itest.eq.0) close(lfntst) 2934c 2935 return 2936 end 2937 subroutine md_fd(isg,xs,fst,fs,iwg,xw,fwt,fw) 2938c 2939 implicit none 2940c 2941#include "md_common.fh" 2942#include "mafdecls.fh" 2943#include "msgids.fh" 2944c 2945 integer isg(msa),iwg(mwm) 2946 real*8 xs(msa,3),fst(msa,3),fs(msa,3) 2947 real*8 xw(mwm,3,mwa),fwt(mwm,3,mwa),fw(mwm,3,mwa) 2948c 2949 integer i,j,k,lsa 2950 real*8 etott,xsk(3),ft(3,3),dx,xdv 2951c 2952 if(me.eq.0) write(lfnout,1000) 2953 1000 format(//,' FINITE DIFFERENCE SOLUTE FORCES',//, 2954 + ' Atom', 2955 + ' Analytic forces ', 2956 + ' Finite difference forces ', 2957 + ' Deviation ',/, 2958 + ' ', 2959 + ' fx fy fz', 2960 + ' fx fy fz', 2961 + ' dfx dfy dfz',/) 2962c 2963 lpair=.false. 2964 lload=.false. 2965c xdv=0.000001 2966 xdv=dx0sd 2967c 2968 do 1 i=1,nsaloc 2969 fst(i,1)=fs(i,1) 2970 fst(i,2)=fs(i,2) 2971 fst(i,3)=fs(i,3) 2972 1 continue 2973 do 2 i=1,nwmloc 2974 do 3 j=1,nwa 2975 fwt(i,1,j)=fw(i,1,j) 2976 fwt(i,2,j)=fw(i,2,j) 2977 fwt(i,3,j)=fw(i,3,j) 2978 3 continue 2979 2 continue 2980 etott=etot 2981c 2982 do 4 i=1,nsa 2983 lsa=0 2984 ft(1,1)=zero 2985 ft(2,1)=zero 2986 ft(3,1)=zero 2987 ft(1,2)=zero 2988 ft(2,2)=zero 2989 ft(3,2)=zero 2990 ft(1,3)=zero 2991 ft(2,3)=zero 2992 ft(3,3)=zero 2993 do 5 j=1,nsaloc 2994 if(isg(j).eq.i) then 2995 lsa=j 2996 xsk(1)=xs(j,1) 2997 xsk(2)=xs(j,2) 2998 xsk(3)=xs(j,3) 2999 ft(1,1)=fst(j,1) 3000 ft(2,1)=fst(j,2) 3001 ft(3,1)=fst(j,3) 3002 goto 6 3003 endif 3004 5 continue 3005 6 continue 3006c 3007 do 7 j=2,3 3008 do 8 k=1,3 3009c 3010 dx=-xdv 3011 if(j.eq.3) dx=xdv 3012c 3013 if(lsa.gt.0) xs(lsa,k)=xs(lsa,k)+dx 3014c 3015c atomic forces and potential energies 3016c 3017 call md_finit(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 3018 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs), 3019 + dbl_mb(i_xsm),dbl_mb(i_xsmp)) 3020 call md_forces(int_mb(i_iw),int_mb(i_is),dbl_mb(i_xw), 3021 + dbl_mb(i_xwm),dbl_mb(i_xs),dbl_mb(i_fw),dbl_mb(i_fs)) 3022c 3023 call prp_proper(0,stime,eww,dbl_mb(i_esw), 3024 + dbl_mb(i_ess),dbl_mb(i_fss),dbl_mb(i_esk),epme,uqmd,uqmmm, 3025 + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot, 3026 + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm)) 3027c 3028c 3029 if(lsa.gt.0) then 3030 xs(lsa,k)=xsk(k) 3031 ft(k,j)=etott-etot 3032 endif 3033c 3034 8 continue 3035 7 continue 3036c 3037 if(np.gt.1) call ga_dgop(mrg_d14,ft,6,'+') 3038c 3039 if(me.eq.0) write(lfnout,1001) i,(ft(j,1),j=1,3), 3040 + ((ft(j,3)-ft(j,2))/(two*xdv),j=1,3), 3041 + (ft(j,1)-(ft(j,3)-ft(j,2))/(two*xdv),j=1,3) 3042 1001 format(i7,5x,3f12.3,3x,3f12.3,3x,3E12.3) 3043c 3044 4 continue 3045c 3046 return 3047 end 3048 subroutine md_membrane_init(ismol,mm,xs,xsm,fm) 3049c 3050 implicit none 3051c 3052#include "md_common.fh" 3053#include "mafdecls.fh" 3054#include "msgids.fh" 3055c 3056 integer ismol(msa),mm(msa,2) 3057 real*8 xs(msa,3),xsm(msm,3),fm(msm,7) 3058c 3059 integer i,j,k,numg 3060 real*8 dx,dk 3061c 3062 numg=0 3063c 3064 do 1 i=1,msm 3065 mm(i,1)=0 3066 1 continue 3067 do 2 i=1,nsaloc 3068 mm(ismol(i),1)=mm(ismol(i),1)+1 3069 mm(i,2)=ismol(i) 3070 2 continue 3071 if(np.gt.1) call ga_dgop(mrg_d49,mm,msm,'+') 3072 do 3 i=1,nsaloc 3073 if(mm(ismol(i),1).eq.1) then 3074 k=0 3075 dk=zero 3076 do 4 j=1,nsm 3077 if(ismol(i).ne.j) then 3078 dx=(xs(i,1)-xsm(j,1))**2+(xs(i,2)-xsm(j,2))**2+ 3079 + (xs(i,3)-xsm(j,3))**2 3080 if(k.eq.0.or.dx.lt.dk) then 3081 dk=dx 3082 k=j 3083 endif 3084 endif 3085 4 continue 3086 mm(i,2)=k 3087 numg=numg+1 3088 endif 3089 3 continue 3090 if(numg.gt.0.and.me.eq.0) then 3091 write(*,2000) numg 3092 2000 format(' Regrouping of',i5,' atoms',/) 3093 endif 3094c 3095 return 3096 end 3097 subroutine md_membrane_forces(mm,fm,xs,xsm,fs,ws) 3098c 3099 implicit none 3100c 3101#include "md_common.fh" 3102#include "mafdecls.fh" 3103#include "msgids.fh" 3104c 3105 integer mm(msa,2) 3106 real*8 fm(msm,7),xs(msa,3),xsm(msm,3),fs(msa,3),ws(msa) 3107c 3108 integer i,j 3109 real*8 factor 3110c 3111 do 1 i=1,msm 3112 fm(i,1)=zero 3113 fm(i,2)=zero 3114 fm(i,3)=zero 3115 fm(i,4)=zero 3116 fm(i,5)=zero 3117 fm(i,6)=zero 3118 fm(i,7)=zero 3119 1 continue 3120c 3121 do 2 i=1,nsaloc 3122 factor=one/ws(i) 3123 fm(mm(i,2),1)=fm(mm(i,2),1)+factor*fs(i,1) 3124 fm(mm(i,2),2)=fm(mm(i,2),2)+factor*fs(i,2) 3125 fm(mm(i,2),3)=fm(mm(i,2),3)+factor*fs(i,3) 3126 fm(mm(i,2),4)=fm(mm(i,2),4)+factor* 3127 + ((xs(i,1)-xsm(mm(i,2),1))*fs(i,2)- 3128 + (xs(i,2)-xsm(mm(i,2),2))*fs(i,1)) 3129 2 continue 3130 if(np.gt.1) call ga_dgop(mrg_d50,fm,4*msm,'+') 3131c 3132 do 3 i=1,nsm 3133 fm(i,4)=fm(i,3) 3134 if(me.eq.0) write(lfnout,1000) i,(fm(i,j)/dble(mm(i,1)),j=1,3) 3135 1000 format(i5,3f12.3) 3136 3 continue 3137c 3138c molecular rotations only 3139c 3140 if(imembr.eq.2) then 3141 do 4 i=1,nsaloc 3142 fs(i,1)=fs(i,1)-ws(i)*fm(mm(i,2),1)/dble(mm(i,1)) 3143 fs(i,2)=fs(i,2)-ws(i)*fm(mm(i,2),2)/dble(mm(i,1)) 3144 fs(i,3)=fs(i,3)-ws(i)*fm(mm(i,2),3)/dble(mm(i,1)) 3145 4 continue 3146 endif 3147c 3148 return 3149 end 3150