1 subroutine md_start 2c 3c $Id$ 4c 5 implicit none 6c 7#include "md_common.fh" 8#include "mafdecls.fh" 9#include "rtdb.fh" 10c 11 character*5 gid 12 integer iun,ibl 13 logical lexist 14c 15 if(npg.gt.1) call md_partition 16c 17 if(me.eq.0) then 18 if(npg.gt.1) then 19 write(gid,3300) ipg 20 3300 format(i5.5) 21 ibl=index(filnam,' ')-1 22 if(ibl.lt.1) then 23 filnam='nwmd' 24 ibl=4 25 endif 26 iun=index(filnam,'_')-1 27 if(iun.lt.0) iun=ibl 28 if(iun.eq.0) then 29 ibl=index(filnam,' ')-1 30 iun=5 31 endif 32 filtop=filnam(1:iun)//'.top' 33 filrst=filnam(1:ibl)//'.rst' 34 filout=filnam(1:ibl)//gid//'.out' 35 filtrj=filnam(1:ibl)//gid//'.trj' 36 filprp=filnam(1:ibl)//gid//'.prp' 37 filmro=filnam(1:ibl)//gid//'.mro' 38 filgib=filnam(1:ibl)//gid//'.gib' 39 if(npg.gt.1) then 40 filrst=filnam(1:ibl)//gid//'.rst' 41 inquire(file=filrst(1:index(filrst,' ')-1),exist=lexist) 42 if(.not.lexist) filrst=filnam(1:ibl)//'.rst' 43 endif 44 endif 45 endif 46c 47c Open debug file if requested 48c ---------------------------- 49c 50 if(idebug.gt.0) then 51 lfndbg=18 52 write(fildbg,'(a,i5.5,a)') 'nwchem_',me,'.dbg' 53 open(unit=lfndbg,file=fildbg,form='formatted',status='unknown') 54 endif 55c 56c open recording file 57c 58 if(me.eq.0) then 59 call md_fopen(.false.) 60c 61 if(ntype.eq.3) then 62 open(unit=lfngib,file=filgib(1:index(filgib,' ')-1), 63 + form='formatted',status='unknown') 64 endif 65 if(ntype.eq.0.and.nftri.gt.0) then 66 open(unit=lfntri,file=filtri(1:index(filtri,' ')-1), 67 + form='formatted',status='unknown') 68 endif 69c 70 if(itest.gt.0) then 71 open(unit=lfntst,file=filtst(1:index(filtst,' ')-1), 72 + form='formatted',status='unknown') 73 endif 74 endif 75c 76 if(lqmd) call qmd_start() 77c 78c print input information 79c ----------------------- 80c 81 call md_print() 82c 83c start the spacial decomposition API 84c ----------------------------------- 85c 86 if(.not.lqmd) then 87 call sp_start(lfnout,lfntop,filtop,lfnrst,filrst,lfnsyn,filsyn, 88 + nfsync,rshort,rlong,zero,rsgm, 89 + npx,npy,npz,nbx,nby,nbz, 90 + npbtyp,nbxtyp,box,vlat,lpbc, 91 + nwm,mwm,nwa,mwa,nsf,msf,nsm,msm,nsa,msa, 92 + loadb,lbpair,factld,ipolt.ne.0,.false.,temp,tempw,temps,lqmd, 93 + iguide,lfndbg,idebug,projct,mbbreq,nserie,isload,ireset,icntrl, 94 + nseq,i_lseq,ndums,nbget,nprec,madbox) 95 endif 96c 97c start the analysis API 98c ---------------------- 99c 100c if(nfanal.gt.0) call ana_init(nsa,msa,.false.) 101c 102c start the classical forces API 103c ------------------------------ 104c 105 call cf_start(irtdb,lqmd,lfnout,lfntop,filtop,ndistr,npmf,npmfi, 106 + nwm,mwm,nwa,mwa,nsf,msf,nsm,msm,nsa,msa, 107 + mdalgo,npbtyp,nbxtyp,rshort,rlong,rqmmm,box, 108 + ipme,morder,ngx,ngy,ngz,nodpme,pmetol, 109 + tstep,tlwsha,mshitw,mshits,noshak, 110 + iqmmm,ipolt,itscal,ipscal,ipopt,tmpext,prsext, 111 + tmprlx,tmsrlx,prsrlx,scaleq,facpmf, 112 + 0,temp,tempw,temps,compr,ntype,iset,isetp1,isetp2, 113 + issscl,delta,nfanal,lpbc,npgdec,xfield,xfvect,xffreq, 114 + npener,icntrl,nbias,mropt,includ,ltwin, 115 + nseq,i_lseq,nfhop,rhop,thop,ndums,ipbtyp,lfnhop,iradgy, 116 + 0,nprec) 117c + nbget) 118c 119 if(mlambd.gt.0.and.ilambd.gt.0.and.ilambd.le.mlambd) 120 + call cf_lambda(lamtyp,irun,maxlam,elam,lfnout,lfnpmf, 121 + rlambd,dlambd,filnam) 122c 123c print topology information 124c -------------------------- 125c 126 call cf_print_top(lfnout,npatom,nptopw,nptops) 127c 128 msm=max(1,nsm) 129 msf=max(1,nsf) 130 mst=max(msm,nseq) 131c 132c start the property API 133c ---------------------- 134c 135 call prp_start(nserie,ntype,nftri,lfnrst,filrst,lfnout,lfnprp, 136 + lfngib,nfoutp,nfstat,nfprop,iprop, 137 + .true.,.true.,ltwin,ipolt.ne.0,ipme.ne.0, 138 + npstep.ne.0,npener.ne.0,npstat, 139 + nwm,msf,nsf,mpe,mdacq,mrun,iset,isetp1,isetp2,tstep,msm,nsm, 140 + nsa,ddacq,edacq,iprof,npmf,npener,ndistr,lpbc,nbias,nodpme,npmfi, 141 + iguide.ne.0,iqmmm.ne.0,lqmd,iradgy,idifco,nbget,ipg,npg) 142c 143c allocate memory for coordinates, velocities, etc. 144c ------------------------------------------------- 145c 146 if(.not.ma_push_get(mt_int,mwm*miw2,'iw',l_iw,i_iw)) 147 + call md_abort('Failed to allocate memory for iw',0) 148 if(.not.ma_push_get(mt_int,msa*mis2,'is',l_is,i_is)) 149 + call md_abort('Failed to allocate memory for is',0) 150 if(.not.ma_push_get(mt_int,mwm,'iwz',l_iwz,i_iwz)) 151 + call md_abort('Failed to allocate memory for iw',0) 152 if(.not.ma_push_get(mt_int,msa,'isz',l_isz,i_isz)) 153 + call md_abort('Failed to allocate memory for is',0) 154 if(.not.ma_push_get(mt_dbl,3*mwa*mwm,'xw',l_xw,i_xw)) 155 + call md_abort('Failed to allocate memory for xw',0) 156 if(.not.ma_push_get(mt_dbl,3*mwm,'xwm',l_xwm,i_xwm)) 157 + call md_abort('Failed to allocate memory for xwm',0) 158 if(.not.ma_push_get(mt_dbl,mwm,'rtos',l_rtos,i_rtos)) 159 + call md_abort('Failed to allocate memory for rtos',0) 160 if(.not.ma_push_get(mt_dbl,3*msa,'xs',l_xs,i_xs)) 161 + call md_abort('Failed to allocate memory for xs',0) 162 if(.not.ma_push_get(mt_dbl,3*mwa*mwm,'yw',l_yw,i_yw)) 163 + call md_abort('Failed to allocate memory for yw',0) 164 if(.not.ma_push_get(mt_dbl,3*msa,'ys',l_ys,i_ys)) 165 + call md_abort('Failed to allocate memory for ys',0) 166 if(.not.ma_push_get(mt_dbl,3*mwa*mwm,'vw',l_vw,i_vw)) 167 + call md_abort('Failed to allocate memory for vw',0) 168 if(.not.ma_push_get(mt_dbl,3*msa,'vs',l_vs,i_vs)) 169 + call md_abort('Failed to allocate memory for vs',0) 170 if(.not.ma_push_get(mt_dbl,3*mwa*mwm,'vwt',l_vwt,i_vwt)) 171 + call md_abort('Failed to allocate memory for vwt',0) 172 if(.not.ma_push_get(mt_dbl,3*msa,'vst',l_vst,i_vst)) 173 + call md_abort('Failed to allocate memory for vst',0) 174 if(.not.ma_push_get(mt_dbl,6*mwa*mwm,'fw',l_fw,i_fw)) 175 + call md_abort('Failed to allocate memory for fw',0) 176 if(.not.ma_push_get(mt_dbl,6*msa,'fs',l_fs,i_fs)) 177 + call md_abort('Failed to allocate memory for fs',0) 178 if(.not.ma_push_get(mt_dbl,3*mwm,'xwcr',l_xwcr,i_xwcr)) 179 + call md_abort('Failed to allocate memory for xwcr',0) 180 if(.not.ma_push_get(mt_dbl,3*msm,'xsm',l_xsm,i_xsm)) 181 + call md_abort('Failed to allocate memory for xsm',0) 182 if(.not.ma_push_get(mt_dbl,4*mst,'tsm',l_tsm,i_tsm)) 183 + call md_abort('Failed to allocate memory for tsm',0) 184 if(.not.ma_push_get(mt_dbl,3*msm,'xsm',l_xsmp,i_xsmp)) 185 + call md_abort('Failed to allocate memory for xsmp',0) 186 if(.not.ma_push_get(mt_dbl,8*msm,'gsm',l_gsm,i_gsm)) 187 + call md_abort('Failed to allocate memory for gsm',0) 188 if(.not.ma_push_get(mt_dbl,3*msm,'xscr',l_xscr,i_xscr)) 189 + call md_abort('Failed to allocate memory for xscr',0) 190 if(.not.ma_push_get(mt_dbl,msm,'dsr',l_dsr,i_dsr)) 191 + call md_abort('Failed to allocate memory for dsr',0) 192 if(.not.ma_push_get(mt_dbl,18*msm,'zs',l_zs,i_zs)) 193 + call md_abort('Failed to allocate memory for zs',0) 194 if(.not.ma_push_get(mt_dbl,2*mpe*msf,'esw',l_esw,i_esw)) 195 + call md_abort('Failed to allocate memory for esw',0) 196 if(.not.ma_push_get(mt_dbl,2*mpe*msf*msf,'ess',l_ess,i_ess)) 197 + call md_abort('Failed to allocate memory for ess',0) 198 if(.not.ma_push_get(mt_dbl,6*msf*msf,'fss',l_fss,i_fss)) 199 + call md_abort('Failed to allocate memory for fss',0) 200 if(.not.ma_push_get(mt_dbl,msf,'esk',l_esk,i_esk)) 201 + call md_abort('Failed to allocate memory for esk',0) 202 if(.not.ma_push_get(mt_dbl,mwa+msa,'wws',l_wws,i_wws)) 203 + call md_abort('Failed to allocate memory for wws',me) 204 if(npener.eq.0) then 205 if(.not.ma_push_get(mt_dbl,1,'esa',l_esa,i_esa)) 206 + call md_abort('Failed to allocate memory for esa',0) 207 else 208 if(.not.ma_push_get(mt_dbl,2*nsa,'esa',l_esa,i_esa)) 209 + call md_abort('Failed to allocate memory for esa',0) 210 endif 211 if(ipolt.gt.0) then 212 if(.not.ma_push_get(mt_dbl,6*mwa*mwm,'pw',l_pw,i_pw)) 213 + call md_abort('Failed to allocate memory for pw',0) 214 if(.not.ma_push_get(mt_dbl,6*msa,'ps',l_ps,i_ps)) 215 + call md_abort('Failed to allocate memory for ps',0) 216 if(lpert2.or.lpert3) then 217 if(.not.ma_push_get(mt_dbl,12*mwa*mwm,'pwp',l_pwp,i_pwp)) 218 + call md_abort('Failed to allocate memory for pwp',0) 219 if(.not.ma_push_get(mt_dbl,12*msa,'psp',l_psp,i_psp)) 220 + call md_abort('Failed to allocate memory for psp',0) 221 endif 222 endif 223 if(iguide.gt.0) then 224 if(.not.ma_push_get(mt_dbl,3*mwa*mwm,'gw',l_gw,i_gw)) 225 + call md_abort('Failed to allocate memory for gw',0) 226 if(.not.ma_push_get(mt_dbl,3*msa,'gs',l_gs,i_gs)) 227 + call md_abort('Failed to allocate memory for gs',0) 228 endif 229 if(icmopt.gt.0) then 230 if(.not.ma_push_get(mt_dbl,5*msm,'fcm',l_fcm,i_fcm)) 231 + call md_abort('Failed to allocate memory for fcm',0) 232 endif 233c 234 if(imembr.gt.0) then 235 if(.not.ma_push_get(mt_int,2*msa,'mm',l_mm,i_mm)) 236 + call md_abort('Failed to allocate memory for mm',me) 237 if(.not.ma_push_get(mt_dbl,7*msm,'fm',l_fm,i_fm)) 238 + call md_abort('Failed to allocate memory for fm',me) 239 endif 240c 241c retrieve current coordinates for this node 242c ------------------------------------------ 243c 244 if(lqmd) then 245 call qmd_setup(nwmloc, 246 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),dbl_mb(i_gs),nsaloc) 247 else 248 call sp_setup(me,int_mb(i_iw),dbl_mb(i_xw),dbl_mb(i_xwcr), 249 + dbl_mb(i_vw),dbl_mb(i_gw),nwmloc, 250 + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_xscr),dbl_mb(i_vs), 251 + dbl_mb(i_gs),nsaloc,lpack) 252 endif 253 if(.not.lqmd) 254 + call sp_update_i(nsaloc,int_mb(i_is),nwmloc,int_mb(i_iw)) 255c 256c initialize packing 257c ------------------ 258c 259 if(lpack) 260 + call sp_pack_init(int_mb(i_is),nsaloc,int_mb(i_iw),nwmloc) 261c 262c spacial decomposition 263c --------------------- 264c 265 if(.not.lqmd) then 266 call sp_travel(box,dbl_mb(i_xw),dbl_mb(i_vw),dbl_mb(i_xwcr), 267 + dbl_mb(i_gw),int_mb(i_iw),nwmloc,dbl_mb(i_xs),dbl_mb(i_vs), 268 + dbl_mb(i_gs),int_mb(i_is),nsaloc) 269 endif 270c 271c calculate mass factors 272c ---------------------- 273c 274 call cf_weight(nwmloc,nsaloc,int_mb(i_is+(lsatt-1)*msa), 275 + int_mb(i_is+(lsmol-1)*msa),int_mb(i_is+(lshop-1)*msa),wbox) 276c 277c calculate centers of mass 278c ------------------------- 279c 280 call cf_cenmas(nwmloc,dbl_mb(i_xw),dbl_mb(i_xwm),nsaloc, 281 + int_mb(i_is+(lsatt-1)*msa),int_mb(i_is+(lsmol-1)*msa), 282 + dbl_mb(i_xs),dbl_mb(i_xsm),dbl_mb(i_gsm)) 283c 284c fix 285c ------------ 286c 287 if(me.eq.0.and.numfix.gt.0) then 288 open(unit=lfncmd,file=filcmd(1:index(filcmd,' ')-1), 289 + form='formatted',status='old') 290 rewind(lfncmd) 291 endif 292 call cf_fix(lfnout,lfncmd,numfix,int_mb(i_iw+(lwgmn-1)*mwm), 293 + int_mb(i_iw+(lwdyn-1)*mwm),nwmloc, 294 + int_mb(i_is+(lsgan-1)*msa),int_mb(i_is+(lsatt-1)*msa), 295 + int_mb(i_is+(lsdyn-1)*msa),int_mb(i_is+(lssgm-1)*msa),nsaloc, 296 + dbl_mb(i_xwm),dbl_mb(i_xs)) 297 if(me.eq.0.and.numfix.gt.0) then 298 close(unit=lfncmd) 299 endif 300 if(.not.lqmd) 301 + call sp_update_i(nsaloc,int_mb(i_is),nwmloc,int_mb(i_iw)) 302c 303c print decomposition information 304c ------------------------------- 305c 306 if(.not.lqmd) call sp_print() 307c 308 call prp_setup(wbox) 309c 310c write file headers 311c ------------------ 312c 313 if(me.eq.0.and.ntype.ne.3) then 314 if(nfcoor.gt.0.or.nfscoo.gt.0.or.nfvelo.gt.0.or.nfsvel.gt.0) 315 + call cf_trjhdr(lfntrj) 316 endif 317c 318 if(npg.gt.1.and.me.eq.0) then 319 root=filnam(1:ibl)//gid 320 filrst=filnam(1:ibl)//gid//'.rst' 321 rfile=filrst 322 endif 323c 324 if(.not.ma_verify_allocator_stuff()) 325 + call md_abort('ma problems at end of md_start',me) 326 return 327 end 328 subroutine md_rdinp() 329c 330 implicit none 331c 332#include "md_common.fh" 333#include "global.fh" 334#include "msgids.fh" 335#include "rtdb.fh" 336#include "mafdecls.fh" 337#include "nwmd_constants.fh" 338#if defined(NEED_LOC) 339 external loc 340 integer loc 341#endif 342c 343 character*20 operat 344 character*32 theory 345 logical lstate,lstart,lrstrt,lcont 346 integer igt,igr,nbytes,ibl,iun 347 integer niperw,nbperw,istart 348 character*5 gid 349c 350c set process id and number of processes 351c -------------------------------------- 352c 353 me=ga_nodeid() 354 np=ga_nnodes() 355c 356c set logical file numbers 357c 358 lfninp=10 359 lfnout=11 360 lfntop=12 361 lfnrst=13 362 lfntrj=14 363 lfnprp=15 364 lfngib=16 365 lfnqrs=17 366 lfndbg=18 367 lfnmri=19 368 lfnmro=20 369 lfntim=21 370 lfnsyn=22 371c 372 lfnmrr=23 373 lfnarg=24 374 lfnsum=25 375 lfnlog=26 376 lfnind=27 377 lfnsin=28 378 lfnrdi=29 379 lfnfld=30 380 lfnsfl=31 381 lfntst=32 382 lfnrdf=33 383 lfnppd=34 384 lfndip=35 385 lfndef=36 386 lfnhis=37 387 lfnhbo=38 388 lfnkrk=39 389 lfnqsc=40 390 lfnacf=41 391 lfncnv=42 392 lfnfet=43 393 lfnuse=44 394 lfnmsg=45 395 lfnday=46 396 lfncmd=47 397 lfnpmf=48 398 lfntri=49 399 lfnhop=50 400c 401c set defaults 402c 403 titinp(1)='NWChem:MD input' 404 titinp(2)=' ' 405 titinp(3)=' ' 406c 407 call swatch(datinp,timinp) 408c 409c print flags 410c 411 nptopw=0 412 nptops=0 413 nptopt=0 414 npstep=0 415 npstat=0 416 npener=0 417 npforc=0 418 nptmng=0 419 npmemo=0 420c 421c default task is single point energy set 1 422c 423 ntype=0 424 mdtype=1 425 iset=1 426 mdordr=0 427 nserie=0 428c 429c polarization defaults 430c 431 ipolt=0 432 mpolit=1 433 ptol=1.0d-5 434c 435c default cutoff radii 436c 437 iswtch=0 438 rshort=0.9d0 439 rlong=0.9d0 440c 441c md defaults 442c 443 stime=zero 444 mdstep=0 445 tstep=0.001d0 446 kequi=0 447 mequi=0 448 kdacq=0 449 mdacq=100 450 npgdrv=0 451 npgdec=0 452 iffdrv=0 453 isaltb=0 454c 455c mcti defaults 456c 457 maxlam=21 458 elam=one 459 ddacq=5.0d0 460 edacq=5.0d0 461 fdacq=0.75d0 462 macfl=1 463 ixcl=0 464 iapprx=0 465 weight=-2.5d-2 466 facapp=two 467c 468 dgscl=one 469 dgref=5.0d-4 470 ddgscl=5.0d-4 471c 472 ssshft=0.075d0 473c 474c multiple run defaults 475c 476 mrun=0 477 mropt=0 478 ldacq=0 479 mdopt=0 480 msplit=0 481c 482c em defaults 483c 484 mintyp=0 485 nem=0 486 nemcrt=3 487 nfqrs=25 488c 489 msdit=100 490 dx0sd=0.01d0 491 dxmsd=0.05d0 492 dxsdmx=0.00001d0 493c 494 mcgit=0 495 ncgcy=0 496 dx0cg=0.01d0 497 dxcgmx=0.00001d0 498c 499c shake defaults 500c 501 mshitw=100 502 mshits=100 503 tlwsha=0.001d0 504 tlssha=0.001d0 505 ifss=0 506c 507c pairlist defaults 508c 509 lwtype=0 510c 511c frequency defaults 512c 513 nfoutp=-1 514 nfstat=-1 515 nfrest=0 516 nffree=-1 517c 518 nfcoor=0 519 nfscoo=0 520 nfprop=0 521 nfvelo=0 522 nfsvel=0 523 nfforc=0 524 nfsfor=0 525 nfpold=0 526 nrwrec=0 527 nfindu=0 528 nfsind=0 529c 530 nfcntr=0 531 nfslow=0 532 nfshkw=0 533 nfshks=0 534c 535 nfrdf=0 536 nfdip=0 537 nfkirk=0 538 nkirk=0 539 nfhbo=0 540 drdf=0.0d0 541c 542 lenhis=0 543 lendis=0 544c 545 nfsync=0 546 iprop=0 547c 548 lserver=.false. 549c 550c distributions 551c 552 ngl=0 553 rrdf=zero 554 ngr=0 555 ngc=0 556 ngrww=0 557 ngrsw=0 558 ngrss=0 559c 560 ndip=0 561 rdip=zero 562 rkirk=zero 563c 564 numhis=0 565 numdis=0 566 lnghis=0 567 lngdis=0 568c 569 nfacfa=0 570 nfauto=0 571 nfconv=0 572c 573 nfdip=0 574 nfkirk=0 575c 576 issscl=0 577c 578 nfnewf=0 579 ibatch=0 580c 581c constant p 582c 583 ipscal=0 584 ipopt=0 585 prsext=1.025d5 586 prsrlx=0.5d0 587 compr=4.53d-10 588c 589c constant t 590c 591 itscal=0 592 tmpext=298.15d0 593 tmpext1=298.15d0 594 tmpext2=298.15d0 595 tmprlx=0.4d0 596 tmsrlx=0.4d0 597 tann1=0.0d00 598 tann2=1.0d02 599c 600c velocity reassignment 601c 602 ivreas=0 603 ivopt=0 604 tgauss=298.15d0 605 iseed=12345 606c 607c reaction field 608c 609 ireact=0 610 dielec=one 611c 612c external electric field 613c 614 ifield=0 615 field=zero 616c 617c solute centering 618c 619 nfcntr=0 620c 621c print default 622c 623 nprint=3 624c 625c recording defaults 626c 627 ibinar=0 628c 629c load balancing 630c 631 factld=zero 632c 633c other defaults (some are inactive) 634c 635 icbw=1 636 icbs=1 637 irr=0 638 sil=0.001d0 639 idipol=0 640 iwarn=0 641 irlow=0 642 idebug=0 643 iumbr=0 644 igmass=0 645 lowmcr=0 646 noone=0 647 npolit=0 648c 649 nwarn=0 650 nform=0 651 nopack=0 652c 653 urlow=tiny 654c 655 ignore=0 656 mdo=6 657 mwork=-1 658 hbdis1=pthree 659 hbdis2=pfour 660 hbang1=two*atan(one) 661 hbang2=two*atan(one) 662 verinp=zero 663c 664 uqmd=zero 665c 666 iformt=1 667c 668c parallel execution defaults 669c 670c many of the above defaults are superceded by defaults in the rtdb 671c 672c general initializations 673c 674c nbperw : number of bytes per word 675c niperw : number of integers per word 676c 677 nbperw=ma_sizeof(mt_dbl,1,mt_byte) 678 niperw=nbperw/ma_sizeof(mt_int,1,mt_byte) 679c 680c get input variables on process 0 only 681c 682 if(me.eq.0) then 683c 684c sequential access to rtdb 685c 686 lstate=rtdb_parallel(.false.) 687c 688c retrieve project name 689c 690 if(.not.rtdb_cget(irtdb,'md:project',1,filnam)) then 691 if(.not.rtdb_cget(irtdb,'prep:sysnam',1,filnam)) then 692 if(.not.rtdb_cget(irtdb,'file_prefix',1,filnam)) 693 + call md_abort('md_rdname: error rtdb_cget project',0) 694 endif 695 endif 696 projct=filnam 697c 698 ibl=index(filnam,' ')-1 699 if(ibl.lt.1) then 700 filnam='nwmd' 701 ibl=4 702 endif 703 iun=index(filnam,'_')-1 704 if(iun.lt.0) iun=ibl 705 if(iun.eq.0) then 706 ibl=index(filnam,' ')-1 707 iun=5 708 endif 709c 710c set file names 711c 712 root=filnam(1:ibl) 713 filday=filnam(1:ibl)//'.day' 714 filinp=filnam(1:ibl)//'.inp' 715 filrst=filnam(1:ibl)//'.rst' 716 filtop=filnam(1:iun)//'.top' 717 filtrj=filnam(1:ibl)//'.trj' 718 filprp=filnam(1:ibl)//'.prp' 719 filgib=filnam(1:ibl)//'.gib' 720 filppd=filnam(1:ibl)//'.ppd' 721 filhis=filnam(1:ibl)//'.his' 722 filhbo=filnam(1:ibl)//'.hbo' 723 filmro=filnam(1:ibl)//'.mro' 724 filmrr=filnam(1:ibl)//'.mrr' 725 filmri=filnam(1:ibl)//'.mri' 726 filout=filnam(1:ibl)//'.out' 727 filmem=filnam(1:ibl)//'.mem' 728 fildef=filnam(1:ibl)//'.def' 729 filqsc=filnam(1:ibl)//'.qsc' 730 filsur=filnam(1:ibl)//'.sur' 731 filqrs=filnam(1:ibl)//'.qrs' 732 filrdf=filnam(1:ibl)//'.rdf' 733 filddf=filnam(1:ibl)//'.ddf' 734 filkrk=filnam(1:ibl)//'.krk' 735 fildbg=filnam(1:ibl)//'.dbg' 736 filacf=filnam(1:ibl)//'.acf' 737 filcnv=filnam(1:ibl)//'.cnv' 738 filfet=filnam(1:ibl)//'.fet' 739 filind=filnam(1:ibl)//'.ind' 740 filsin=filnam(1:ibl)//'.sin' 741 filtim=filnam(1:ibl)//'.tim' 742 filsyn=filnam(1:ibl)//'.syn' 743 filrdi=filnam(1:ibl)//'.rdi' 744 filfld=filnam(1:ibl)//'.fld' 745 filsfl=filnam(1:ibl)//'.sfl' 746 filtst=filnam(1:ibl)//'.tst' 747 filcmd=filnam(1:ibl)//'.cmd' 748 filtri=filnam(1:ibl)//'.tri' 749 filpmf=filnam(1:ibl)//'.pmf' 750 filhop=filnam(1:ibl)//'.hop' 751 rfile=filrst 752c 753c input from rtdb 754c 755 if(.not.rtdb_get(irtdb,'md:port',mt_int,1,iport)) 756 + call md_abort('md_rtdbin: rtdb_get failed',0) 757 if(iport.gt.0) then 758 if(.not.rtdb_cget(irtdb,'md:server',1,server)) 759 + call md_abort('md_rtdbin: rtdb_get failed',0) 760 endif 761 if(.not.rtdb_get(irtdb,'md:istart',mt_int,1,istart)) 762 + call md_abort('md_rtdbin: rtdb_get failed',1) 763 if(.not.rtdb_get(irtdb,'md:nbx',mt_int,1,nbx)) 764 + call md_abort('md_rtdbin: rtdb_get failed',1) 765 if(.not.rtdb_get(irtdb,'md:nby',mt_int,1,nby)) 766 + call md_abort('md_rtdbin: rtdb_get failed',2) 767 if(.not.rtdb_get(irtdb,'md:nbz',mt_int,1,nbz)) 768 + call md_abort('md_rtdbin: rtdb_get failed',3) 769 if(.not.rtdb_get(irtdb,'md:npx',mt_int,1,npx)) 770 + call md_abort('md_rtdbin: rtdb_get failed',4) 771 if(.not.rtdb_get(irtdb,'md:npy',mt_int,1,npy)) 772 + call md_abort('md_rtdbin: rtdb_get failed',5) 773 if(.not.rtdb_get(irtdb,'md:npz',mt_int,1,npz)) 774 + call md_abort('md_rtdbin: rtdb_get failed',6) 775 if(.not.rtdb_get(irtdb,'md:npg',mt_int,1,npg)) 776 + call md_abort('md_rtdbin: rtdb_get failed',6) 777 if(.not.rtdb_get(irtdb,'md:mdalgo',mt_int,1,mdalgo)) 778 + call md_abort('md_rtdbin: rtdb_get failed',7) 779 if(.not.rtdb_get(irtdb,'md:iset',mt_int,1,iset)) 780 + call md_abort('md_rtdbin: rtdb_get failed',8) 781 if(.not.rtdb_get(irtdb,'md:isetp1',mt_int,1,isetp1)) 782 + call md_abort('md_rtdbin: rtdb_get failed',9) 783 if(.not.rtdb_get(irtdb,'md:isetp2',mt_int,1,isetp2)) 784 + call md_abort('md_rtdbin: rtdb_get failed',10) 785 if(.not.rtdb_get(irtdb,'md:iforw',mt_int,1,lamtyp)) 786 + call md_abort('md_rtdbin: rtdb_get failed',11) 787 if(.not.rtdb_get(irtdb,'md:msdit',mt_int,1,msdit)) 788 + call md_abort('md_rtdbin: rtdb_get failed',12) 789 if(.not.rtdb_get(irtdb,'md:mcgit',mt_int,1,mcgit)) 790 + call md_abort('md_rtdbin: rtdb_get failed',13) 791 if(.not.rtdb_get(irtdb,'md:icrit',mt_int,1,icrit)) 792 + call md_abort('md_rtdbin: rtdb_get failed',13) 793 if(.not.rtdb_get(irtdb,'md:ncgcy',mt_int,1,ncgcy)) 794 + call md_abort('md_rtdbin: rtdb_get failed',14) 795 if(.not.rtdb_get(irtdb,'md:mrun',mt_int,1,mrun)) 796 + call md_abort('md_rtdbin: rtdb_get failed',15) 797 if(.not.rtdb_get(irtdb,'md:maxlam',mt_int,1,maxlam)) 798 + call md_abort('md_rtdbin: rtdb_get failed',16) 799 if(.not.rtdb_get(irtdb,'md:npgdec',mt_int,1,npgdec)) 800 + call md_abort('md_rtdbin: rtdb_get failed',17) 801 if(.not.rtdb_get(irtdb,'md:issscl',mt_int,1,issscl)) 802 + call md_abort('md_rtdbin: rtdb_get failed',18) 803 if(.not.rtdb_get(irtdb,'md:mequi',mt_int,1,mequi)) 804 + call md_abort('md_rtdbin: rtdb_get failed',19) 805 if(.not.rtdb_get(irtdb,'md:mdacq',mt_int,1,mdacq)) 806 + call md_abort('md_rtdbin: rtdb_get failed',20) 807 if(.not.rtdb_get(irtdb,'md:ldacq',mt_int,1,ldacq)) 808 + call md_abort('md_rtdbin: rtdb_get failed',21) 809 if(.not.rtdb_get(irtdb,'md:iapprx',mt_int,1,iapprx)) 810 + call md_abort('md_rtdbin: rtdb_get failed',22) 811 if(.not.rtdb_get(irtdb,'md:nacfl',mt_int,1,nacfl)) 812 + call md_abort('md_rtdbin: rtdb_get failed',23) 813 if(.not.rtdb_get(irtdb,'md:ipscal',mt_int,1,ipscal)) 814 + call md_abort('md_rtdbin: rtdb_get failed',24) 815 if(.not.rtdb_get(irtdb,'md:ipopt',mt_int,1,ipopt)) 816 + call md_abort('md_rtdbin: rtdb_get failed',24) 817 if(.not.rtdb_get(irtdb,'md:ivopt',mt_int,1,ivopt)) 818 + call md_abort('md_rtdbin: rtdb_get failed',24) 819 if(.not.rtdb_get(irtdb,'md:itscal',mt_int,1,itscal)) 820 + call md_abort('md_rtdbin: rtdb_get failed',25) 821 if(.not.rtdb_get(irtdb,'md:nfgaus',mt_int,1,nfgaus)) 822 + call md_abort('md_rtdbin: rtdb_get failed',26) 823 if(.not.rtdb_get(irtdb,'md:ipolt',mt_int,1,ipolt)) 824 + call md_abort('md_rtdbin: rtdb_get failed',27) 825 if(.not.rtdb_get(irtdb,'md:mpolit',mt_int,1,mpolit)) 826 + call md_abort('md_rtdbin: rtdb_get failed',28) 827 if(.not.rtdb_get(irtdb,'md:mshitw',mt_int,1,mshitw)) 828 + call md_abort('md_rtdbin: rtdb_get failed',29) 829 if(.not.rtdb_get(irtdb,'md:mshits',mt_int,1,mshits)) 830 + call md_abort('md_rtdbin: rtdb_get failed',30) 831 if(.not.rtdb_get(irtdb,'md:itrack',mt_int,1,itrack)) 832 + call md_abort('md_rtdbin: rtdb_get failed',31) 833 if(.not.rtdb_get(irtdb,'md:npstep',mt_int,1,npstep)) 834 + call md_abort('md_rtdbin: rtdb_get failed',32) 835 if(.not.rtdb_get(irtdb,'md:npstat',mt_int,1,npstat)) 836 + call md_abort('md_rtdbin: rtdb_get failed',32) 837 if(.not.rtdb_get(irtdb,'md:npener',mt_int,1,npener)) 838 + call md_abort('md_rtdbin: rtdb_get failed',33) 839 if(.not.rtdb_get(irtdb,'md:npforc',mt_int,1,npforc)) 840 + call md_abort('md_rtdbin: rtdb_get failed',33) 841 if(.not.rtdb_get(irtdb,'md:npdist',mt_int,1,npdist)) 842 + call md_abort('md_rtdbin: rtdb_get failed',34) 843 if(.not.rtdb_get(irtdb,'md:nptmng',mt_int,1,nptmng)) 844 + call md_abort('md_rtdbin: rtdb_get failed',35) 845 if(.not.rtdb_get(irtdb,'md:npatom',mt_int,1,npatom)) 846 + call md_abort('md_rtdbin: rtdb_get failed',36) 847 if(.not.rtdb_get(irtdb,'md:nptopw',mt_int,1,nptopw)) 848 + call md_abort('md_rtdbin: rtdb_get failed',37) 849 if(.not.rtdb_get(irtdb,'md:nptops',mt_int,1,nptops)) 850 + call md_abort('md_rtdbin: rtdb_get failed',38) 851 if(.not.rtdb_get(irtdb,'md:npxpct',mt_int,1,npxpct)) 852 + call md_abort('md_rtdbin: rtdb_get failed',39) 853 if(.not.rtdb_get(irtdb,'md:nfpair',mt_int,1,nfpair)) 854 + call md_abort('md_rtdbin: rtdb_get failed',40) 855 if(.not.rtdb_get(irtdb,'md:nfesp',mt_int,1,nfesp)) 856 + call md_abort('md_rtdbin: rtdb_get failed',40) 857 if(.not.rtdb_get(irtdb,'md:nfrdf',mt_int,1,nfrdf)) 858 + call md_abort('md_rtdbin: rtdb_get failed',41) 859 if(.not.rtdb_get(irtdb,'md:nflong',mt_int,1,nflong)) 860 + call md_abort('md_rtdbin: rtdb_get failed',42) 861 if(.not.rtdb_get(irtdb,'md:nfcntr',mt_int,1,nfcntr)) 862 + call md_abort('md_rtdbin: rtdb_get failed',43) 863 if(.not.rtdb_get(irtdb,'md:icentr',mt_int,1,icentr)) 864 + call md_abort('md_rtdbin: rtdb_get failed',43) 865 if(.not.rtdb_get(irtdb,'md:nfanal',mt_int,1,nfanal)) 866 + call md_abort('md_rtdbin: rtdb_get failed',43) 867 if(.not.rtdb_get(irtdb,'md:nfslow',mt_int,1,nfslow)) 868 + call md_abort('md_rtdbin: rtdb_get failed',44) 869 if(.not.rtdb_get(irtdb,'md:nfoutp',mt_int,1,nfoutp)) 870 + call md_abort('md_rtdbin: rtdb_get failed',45) 871 if(.not.rtdb_get(irtdb,'md:nfstat',mt_int,1,nfstat)) 872 + call md_abort('md_rtdbin: rtdb_get failed',46) 873 if(.not.rtdb_get(irtdb,'md:nfrest',mt_int,1,nfrest)) 874 + call md_abort('md_rtdbin: rtdb_get failed',47) 875 if(.not.rtdb_get(irtdb,'md:keepr',mt_int,1,keepr)) 876 + call md_abort('md_rtdbin: rtdb_get failed',48) 877 if(.not.rtdb_get(irtdb,'md:nfcoor',mt_int,1,nfcoor)) 878 + call md_abort('md_rtdbin: rtdb_get failed',50) 879 if(.not.rtdb_get(irtdb,'md:nfscoo',mt_int,1,nfscoo)) 880 + call md_abort('md_rtdbin: rtdb_get failed',51) 881 if(.not.rtdb_get(irtdb,'md:nfvelo',mt_int,1,nfvelo)) 882 + call md_abort('md_rtdbin: rtdb_get failed',52) 883 if(.not.rtdb_get(irtdb,'md:nfsvel',mt_int,1,nfsvel)) 884 + call md_abort('md_rtdbin: rtdb_get failed',53) 885 if(.not.rtdb_get(irtdb,'md:nfforc',mt_int,1,nfforc)) 886 + call md_abort('md_rtdbin: rtdb_get failed',50) 887 if(.not.rtdb_get(irtdb,'md:nfsfor',mt_int,1,nfsfor)) 888 + call md_abort('md_rtdbin: rtdb_get failed',51) 889 if(.not.rtdb_get(irtdb,'md:nfprop',mt_int,1,nfprop)) 890 + call md_abort('md_rtdbin: rtdb_get failed',54) 891 if(.not.rtdb_get(irtdb,'md:iprop',mt_int,1,iprop)) 892 + call md_abort('md_rtdbin: rtdb_get failed',54) 893 if(.not.rtdb_get(irtdb,'md:nffree',mt_int,1,nffree)) 894 + call md_abort('md_rtdbin: rtdb_get failed',55) 895 if(.not.rtdb_get(irtdb,'md:nfsync',mt_int,1,nfsync)) 896 + call md_abort('md_rtdbin: rtdb_get failed',56) 897 if(.not.rtdb_get(irtdb,'md:nfauto',mt_int,1,nfauto)) 898 + call md_abort('md_rtdbin: rtdb_get failed',57) 899 if(.not.rtdb_get(irtdb,'md:nfconv',mt_int,1,nfconv)) 900 + call md_abort('md_rtdbin: rtdb_get failed',58) 901 if(.not.rtdb_get(irtdb,'md:nffet',mt_int,1,nffet)) 902 + call md_abort('md_rtdbin: rtdb_get failed',59) 903 if(.not.rtdb_get(irtdb,'md:impfr',mt_int,1,impfr)) 904 + call md_abort('md_rtdbin: rtdb_get failed',59) 905 if(.not.rtdb_get(irtdb,'md:impto',mt_int,1,impto)) 906 + call md_abort('md_rtdbin: rtdb_get failed',59) 907 if(.not.rtdb_get(irtdb,'md:nftri',mt_int,1,nftri)) 908 + call md_abort('md_rtdbin: rtdb_get failed',59) 909 if(.not.rtdb_get(irtdb,'md:iformt',mt_int,1,iformt)) 910 + call md_abort('md_rtdbin: rtdb_get failed',60) 911 if(.not.rtdb_get(irtdb,'md:madbox',mt_int,1,madbox)) 912 + call md_abort('md_rtdbin: rtdb_get failed',61) 913 if(.not.rtdb_get(irtdb,'md:loadb',mt_int,1,loadb)) 914 + call md_abort('md_rtdbin: rtdb_get failed',62) 915 if(.not.rtdb_get(irtdb,'md:ireset',mt_int,1,ireset)) 916 + call md_abort('md_rtdbin: rtdb_get failed',62) 917 if(.not.rtdb_get(irtdb,'md:mropt',mt_int,1,mropt)) 918 + call md_abort('md_rtdbin: rtdb_get failed',63) 919 if(.not.rtdb_get(irtdb,'md:idebug',mt_int,1,idebug)) 920 + call md_abort('md_rtdbin: rtdb_get failed',64) 921 if(.not.rtdb_get(irtdb,'md:icntrl',mt_int,1,icntrl)) 922 + call md_abort('md_rtdbin: rtdb_get failed',64) 923 if(.not.rtdb_get(irtdb,'md:ifidi',mt_int,1,ifidi)) 924 + call md_abort('md_rtdbin: rtdb_get failed',64) 925 if(.not.rtdb_get(irtdb,'md:ipbtyp',mt_int,1,ipbtyp)) 926 + call md_abort('md_rtdbin: rtdb_get ipbtyp failed',64) 927 if(.not.rtdb_get(irtdb,'md:ngl',mt_int,1,ngl)) 928 + call md_abort('md_rtdbin: rtdb_get failed',65) 929 if(.not.rtdb_get(irtdb,'md:ifield',mt_int,1,ifield)) 930 + call md_abort('md_rtdbin: rtdb_get failed',66) 931 if(.not.rtdb_get(irtdb,'md:ipme',mt_int,1,ipme)) 932 + call md_abort('md_rtdbin: rtdb_get failed',67) 933 if(.not.rtdb_get(irtdb,'md:ngx',mt_int,1,ngx)) 934 + call md_abort('md_rtdbin: rtdb_get failed',68) 935 if(.not.rtdb_get(irtdb,'md:ngy',mt_int,1,ngy)) 936 + call md_abort('md_rtdbin: rtdb_get failed',69) 937 if(.not.rtdb_get(irtdb,'md:ngz',mt_int,1,ngz)) 938 + call md_abort('md_rtdbin: rtdb_get failed',70) 939 if(.not.rtdb_get(irtdb,'md:numfix',mt_int,1,numfix)) 940 + call md_abort('md_rtdbin: rtdb_get failed',71) 941c if(.not.rtdb_get(irtdb,'md:iunfix',mt_int,1,iunfix)) 942c + call md_abort('md_rtdbin: rtdb_get failed',72) 943c if(.not.rtdb_get(irtdb,'md:lsffix',mt_int,msf,lsffix)) 944c + call md_abort('md_rtdbin: rtdb_get failed',72) 945 if(.not.rtdb_get(irtdb,'md:noshak',mt_int,1,noshak)) 946 + call md_abort('md_rtdbin: rtdb_get failed',73) 947 if(.not.rtdb_get(irtdb,'md:nfefld',mt_int,1,nfefld)) 948 + call md_abort('md_rtdbin: rtdb_get failed',74) 949 if(.not.rtdb_get(irtdb,'md:nfsfld',mt_int,1,nfsfld)) 950 + call md_abort('md_rtdbin: rtdb_get failed',75) 951 if(.not.rtdb_get(irtdb,'md:nscb',mt_int,1,nscb)) 952 + call md_abort('md_rtdbin: rtdb_get failed',76) 953 if(nscb.gt.10) call md_abort('Increase dimension of idscb',0) 954 if(.not.rtdb_get(irtdb,'md:idscb',mt_int,nscb,idscb)) 955 + call md_abort('md_rtdbin: rtdb_get failed',77) 956 if(.not.rtdb_get(irtdb,'md:ireact',mt_int,1,ireact)) 957 + call md_abort('md_rtdbin: rtdb_get failed',78) 958 if(.not.rtdb_get(irtdb,'md:memlim',mt_int,1,memlim)) 959 + call md_abort('md_rtdbin: rtdb_get failed',79) 960 if(.not.rtdb_get(irtdb,'md:morder',mt_int,1,morder)) 961 + call md_abort('md_rtdbin: rtdb_get failed',80) 962 if(.not.rtdb_get(irtdb,'md:isolvo',mt_int,1,isolvo)) 963 + call md_abort('md_rtdbin: rtdb_get failed',81) 964 if(.not.rtdb_get(irtdb,'md:lfout6',mt_int,1,lfout6)) 965 + call md_abort('md_rtdbin: rtdb_get failed',82) 966 if(.not.rtdb_get(irtdb,'md:iprpmf',mt_int,1,iprpmf)) 967 + call md_abort('md_rtdbin: rtdb_get failed',82) 968 if(.not.rtdb_get(irtdb,'md:imfft',mt_int,1,imfft)) 969 + call md_abort('md_rtdbin: rtdb_get failed',83) 970 if(.not.rtdb_get(irtdb,'md:mwmreq',mt_int,1,mwmreq)) 971 + call md_abort('md_rtdbin: rtdb_get failed',84) 972 if(.not.rtdb_get(irtdb,'md:msareq',mt_int,1,msareq)) 973 + call md_abort('md_rtdbin: rtdb_get failed',85) 974 if(.not.rtdb_get(irtdb,'md:mbbreq',mt_int,1,mbbreq)) 975 + call md_abort('md_rtdbin: rtdb_get failed',85) 976 if(.not.rtdb_get(irtdb,'md:itest',mt_int,1,itest)) 977 + call md_abort('md_rtdbin: rtdb_get failed',86) 978 if(.not.rtdb_get(irtdb,'md:nodpme',mt_int,1,nodpme)) 979 + call md_abort('md_rtdbin: rtdb_get failed',87) 980 if(.not.rtdb_get(irtdb,'md:lbpair',mt_int,1,lbpair)) 981 + call md_abort('md_rtdbin: rtdb_get failed',88) 982 if(.not.rtdb_get(irtdb,'md:ndistr',mt_int,1,ndistr)) 983 + call md_abort('md_rtdbin: rtdb_get failed',89) 984 if(.not.rtdb_get(irtdb,'md:npmf',mt_int,1,npmf)) 985 + call md_abort('md_rtdbin: rtdb_get failed',89) 986 if(.not.rtdb_get(irtdb,'md:facpmf',mt_dbl,1,facpmf)) 987 + call md_abort('md_rtdbin: rtdb_get failed',89) 988 if(.not.rtdb_get(irtdb,'md:ndaver',mt_int,1,ndaver)) 989 + call md_abort('md_rtdbin: rtdb_get failed',90) 990 if(.not.rtdb_get(irtdb,'md:idevel',mt_int,1,idevel)) 991 + call md_abort('md_rtdbin: rtdb_get failed',91) 992c if(.not.rtdb_get(irtdb,'md:itime',mt_int,mtimes,itime)) 993c + call md_abort('md_rtdbin: rtdb_get failed',92) 994 if(.not.rtdb_get(irtdb,'md:nftime',mt_int,1,nftime)) 995 + call md_abort('md_rtdbin: rtdb_get failed',92) 996 if(.not.rtdb_get(irtdb,'md:nfdrss',mt_int,1,nfdrss)) 997 + call md_abort('md_rtdbin: rtdb_get failed',93) 998 if(.not.rtdb_get(irtdb,'md:nfload',mt_int,1,nfload)) 999 + call md_abort('md_rtdbin: rtdb_get failed',94) 1000 if(.not.rtdb_get(irtdb,'md:itload',mt_int,1,itload)) 1001 + call md_abort('md_rtdbin: rtdb_get failed',94) 1002 if(.not.rtdb_get(irtdb,'md:ioload',mt_int,1,ioload)) 1003 + call md_abort('md_rtdbin: rtdb_get failed',94) 1004 if(.not.rtdb_get(irtdb,'md:isload',mt_int,1,isload)) 1005 + call md_abort('md_rtdbin: rtdb_get failed',94) 1006 if(.not.rtdb_get(irtdb,'md:ihess',mt_int,1,ihess)) 1007 + call md_abort('md_rtdbin: rtdb_get failed',95) 1008 if(.not.rtdb_get(irtdb,'md:latom',mt_int,1,latom)) 1009 + call md_abort('md_rtdbin: rtdb_get failed',96) 1010 if(.not.rtdb_get(irtdb,'md:icomm',mt_int,1,icomm)) 1011 + call md_abort('md_rtdbin: rtdb_get failed',97) 1012 if(.not.rtdb_get(irtdb,'md:nfhop',mt_int,1,nfhop)) 1013 + call md_abort('md_rtdbin: rtdb_get failed',98) 1014 if(.not.rtdb_get(irtdb,'md:rhop',mt_dbl,1,rhop)) 1015 + call md_abort('md_rtdbin: rtdb_get failed',98) 1016 if(.not.rtdb_get(irtdb,'md:thop',mt_dbl,1,thop)) 1017 + call md_abort('md_rtdbin: rtdb_get failed',98) 1018 if(.not.rtdb_get(irtdb,'md:iguide',mt_int,1,iguide)) 1019 + call md_abort('md_rtdbin: rtdb_get failed',98) 1020 if(.not.rtdb_get(irtdb,'md:icmopt',mt_int,1,icmopt)) 1021 + call md_abort('md_rtdbin: rtdb_get failed',98) 1022 if(.not.rtdb_get(irtdb,'md:imembr',mt_int,1,imembr)) 1023 + call md_abort('md_rtdbin: rtdb_get failed',98) 1024 if(.not.rtdb_get(irtdb,'md:nopack',mt_int,1,nopack)) 1025 + call md_abort('md_rtdbin: rtdb_get failed',99) 1026 if(.not.rtdb_get(irtdb,'md:iprof',mt_int,1,iprof)) 1027 + call md_abort('md_rtdbin: rtdb_get failed',99) 1028 if(.not.rtdb_get(irtdb,'md:ncoll',mt_int,1,ncoll)) 1029 + call md_abort('md_rtdbin: rtdb_get failed',89) 1030 if(.not.rtdb_get(irtdb,'md:ilambd',mt_int,1,ilambd)) 1031 + call md_abort('md_rtdbin: rtdb_get failed',89) 1032 if(.not.rtdb_get(irtdb,'md:mlambd',mt_int,1,mlambd)) 1033 + call md_abort('md_rtdbin: rtdb_get failed',89) 1034 if(.not.rtdb_get(irtdb,'md:includ',mt_int,1,includ)) 1035 + call md_abort('md_rtdbin: rtdb_get failed',89) 1036 if(.not.rtdb_get(irtdb,'md:iradgy',mt_int,1,iradgy)) 1037 + call md_abort('md_rtdbin: rtdb_get failed',89) 1038 if(.not.rtdb_get(irtdb,'md:idifco',mt_int,1,idifco)) 1039 + call md_abort('md_rtdbin: rtdb_get failed',89) 1040 if(.not.rtdb_get(irtdb,'md:nfnewf',mt_int,1,nfnewf)) 1041 + call md_abort('md_rtdbin: rtdb_get failed',89) 1042 if(.not.rtdb_get(irtdb,'md:nbget',mt_int,1,nbget)) 1043 + call md_abort('md_rtdbin: rtdb_get failed',89) 1044 if(.not.rtdb_get(irtdb,'md:nprec',mt_int,1,nprec)) 1045 + call md_abort('md_rtdbin: rtdb_get failed',89) 1046c 1047 if(.not.rtdb_get(irtdb,'md:fguide',mt_dbl,1,fguide)) 1048 + call md_abort('md_rtdbin: rtdb_get failed',100) 1049 if(.not.rtdb_get(irtdb,'md:tguide',mt_dbl,1,tguide)) 1050 + call md_abort('md_rtdbin: rtdb_get failed',100) 1051 if(.not.rtdb_get(irtdb,'md:dx0sd',mt_dbl,1,dx0sd)) 1052 + call md_abort('md_rtdbin: rtdb_get failed',100) 1053 if(.not.rtdb_get(irtdb,'md:dxsdmx',mt_dbl,1,dxsdmx)) 1054 + call md_abort('md_rtdbin: rtdb_get failed',101) 1055 if(.not.rtdb_get(irtdb,'md:dxmsd',mt_dbl,1,dxmsd)) 1056 + call md_abort('md_rtdbin: rtdb_get failed',102) 1057 if(.not.rtdb_get(irtdb,'md:dx0cg',mt_dbl,1,dx0cg)) 1058 + call md_abort('md_rtdbin: rtdb_get failed',103) 1059 if(.not.rtdb_get(irtdb,'md:dxcgmx',mt_dbl,1,dxcgmx)) 1060 + call md_abort('md_rtdbin: rtdb_get failed',104) 1061 if(.not.rtdb_get(irtdb,'md:dxmcg',mt_dbl,1,dxmcg)) 1062 + call md_abort('md_rtdbin: rtdb_get failed',105) 1063 if(.not.rtdb_get(irtdb,'md:edacq',mt_dbl,1,edacq)) 1064 + call md_abort('md_rtdbin: rtdb_get failed',106) 1065 if(.not.rtdb_get(irtdb,'md:ddacq',mt_dbl,1,ddacq)) 1066 + call md_abort('md_rtdbin: rtdb_get failed',107) 1067 if(.not.rtdb_get(irtdb,'md:fdacq',mt_dbl,1,fdacq)) 1068 + call md_abort('md_rtdbin: rtdb_get failed',108) 1069 if(.not.rtdb_get(irtdb,'md:delta',mt_dbl,1,delta)) 1070 + call md_abort('md_rtdbin: rtdb_get failed',109) 1071 if(.not.rtdb_get(irtdb,'md:stime',mt_dbl,1,stime)) 1072 + call md_abort('md_rtdbin: rtdb_get failed',110) 1073 if(.not.rtdb_get(irtdb,'md:tstep',mt_dbl,1,tstep)) 1074 + call md_abort('md_rtdbin: rtdb_get failed',111) 1075 if(.not.rtdb_get(irtdb,'md:prsext',mt_dbl,1,prsext)) 1076 + call md_abort('md_rtdbin: rtdb_get failed',112) 1077 if(.not.rtdb_get(irtdb,'md:prsrlx',mt_dbl,1,prsrlx)) 1078 + call md_abort('md_rtdbin: rtdb_get failed',113) 1079 if(.not.rtdb_get(irtdb,'md:compr',mt_dbl,1,compr)) 1080 + call md_abort('md_rtdbin: rtdb_get failed',114) 1081 if(.not.rtdb_get(irtdb,'md:tmpext',mt_dbl,1,tmpext1)) 1082 + call md_abort('md_rtdbin: rtdb_get failed',115) 1083 tmpext=tmpext1 1084 if(.not.rtdb_get(irtdb,'md:tmpext2',mt_dbl,1,tmpext2)) 1085 + call md_abort('md_rtdbin: rtdb_get failed',115) 1086 if(.not.rtdb_get(irtdb,'md:tmprlx',mt_dbl,1,tmprlx)) 1087 + call md_abort('md_rtdbin: rtdb_get failed',116) 1088 if(.not.rtdb_get(irtdb,'md:tmsrlx',mt_dbl,1,tmsrlx)) 1089 + call md_abort('md_rtdbin: rtdb_get failed',117) 1090 if(.not.rtdb_get(irtdb,'md:tann1',mt_dbl,1,tann1)) 1091 + call md_abort('md_rtdbin: rtdb_get failed',116) 1092 if(.not.rtdb_get(irtdb,'md:tann2',mt_dbl,1,tann2)) 1093 + call md_abort('md_rtdbin: rtdb_get failed',116) 1094 if(.not.rtdb_get(irtdb,'md:tgauss',mt_dbl,1,tgauss)) 1095 + call md_abort('md_rtdbin: rtdb_get failed',118) 1096 if(.not.rtdb_get(irtdb,'md:frgaus',mt_dbl,1,frgaus)) 1097 + call md_abort('md_rtdbin: rtdb_get failed',119) 1098 if(.not.rtdb_get(irtdb,'md:rlong',mt_dbl,1,rlong)) 1099 + call md_abort('md_rtdbin: rtdb_get failed',120) 1100 if(.not.rtdb_get(irtdb,'md:rshort',mt_dbl,1,rshort)) 1101 + call md_abort('md_rtdbin: rtdb_get failed',121) 1102 if(.not.rtdb_get(irtdb,'md:rqmmm',mt_dbl,1,rqmmm)) 1103 + call md_abort('md_rtdbin: rtdb_get failed',122) 1104 if(.not.rtdb_get(irtdb,'md:ptol',mt_dbl,1,ptol)) 1105 + call md_abort('md_rtdbin: rtdb_get failed',123) 1106 if(.not.rtdb_get(irtdb,'md:tlwsha',mt_dbl,1,tlwsha)) 1107 + call md_abort('md_rtdbin: rtdb_get failed',124) 1108 if(.not.rtdb_get(irtdb,'md:tlssha',mt_dbl,1,tlssha)) 1109 + call md_abort('md_rtdbin: rtdb_get failed',125) 1110 if(.not.rtdb_get(irtdb,'md:factld',mt_dbl,1,factld)) 1111 + call md_abort('md_rtdbin: rtdb_get failed',126) 1112 if(.not.rtdb_get(irtdb,'md:rrdf',mt_dbl,1,rrdf)) 1113 + call md_abort('md_rtdbin: rtdb_get failed',127) 1114 if(.not.rtdb_get(irtdb,'md:xfield',mt_dbl,1,xfield)) 1115 + call md_abort('md_rtdbin: rtdb_get failed',128) 1116 if(.not.rtdb_get(irtdb,'md:xffreq',mt_dbl,1,xffreq)) 1117 + call md_abort('md_rtdbin: rtdb_get failed',129) 1118 if(.not.rtdb_get(irtdb,'md:xfvect',mt_dbl,3,xfvect)) 1119 + call md_abort('md_rtdbin: rtdb_get failed',130) 1120 if(.not.rtdb_get(irtdb,'md:weight',mt_dbl,1,weight)) 1121 + call md_abort('md_rtdbin: rtdb_get failed',131) 1122 if(.not.rtdb_get(irtdb,'md:dielec',mt_dbl,1,dielec)) 1123 + call md_abort('md_rtdbin: rtdb_get failed',132) 1124 if(.not.rtdb_get(irtdb,'md:ealpha',mt_dbl,1,ealpha)) 1125 + call md_abort('md_rtdbin: rtdb_get failed',133) 1126 if(.not.rtdb_get(irtdb,'md:rbox',mt_dbl,1,rbox)) 1127 + call md_abort('md_rtdbin: rtdb_get failed',134) 1128 if(.not.rtdb_get(irtdb,'md:rsgm',mt_dbl,1,rsgm)) 1129 + call md_abort('md_rtdbin: rtdb_get failed',134) 1130 if(.not.rtdb_get(irtdb,'md:disrlx',mt_dbl,1,disrlx)) 1131 + call md_abort('md_rtdbin: rtdb_get failed',135) 1132 if(.not.rtdb_get(irtdb,'md:drsscl',mt_dbl,1,drsscl)) 1133 + call md_abort('md_rtdbin: rtdb_get failed',136) 1134 if(.not.rtdb_get(irtdb,'md:pmetol',mt_dbl,1,pmetol)) 1135 + pmetol=1.0d-05 1136 if(.not.rtdb_get(irtdb,'md:fcoll',mt_dbl,1,fcoll)) 1137 + fcoll=1.0d+01 1138 if(.not.rtdb_get(irtdb,'md:scaleq',mt_dbl,1,scaleq)) 1139 + scaleq=-1.0d+00 1140c 1141 if(.not.rtdb_cget(irtdb,'task:theory',1,theory)) 1142 + call md_abort('md_rtdbin: rtdg_cget failed',137) 1143 if(theory.eq.'md') then 1144 iquant=0 1145 iqmmm=0 1146 else 1147 iqmmm=1 1148 if(.not.rtdb_get(irtdb,'task:QMMM',mt_log,1,lqmmm)) iqmmm=0 1149 if(iqmmm.eq.1) then 1150 if(.not.rtdb_get(irtdb,'qmmm:uqmatm',mt_dbl,1,uqmatm)) 1151 + uqmatm=0.0d0 1152 if(.not.rtdb_get(irtdb,'qmmm:linkatm',mt_int,1,linkatm)) 1153 + linkatm=0 1154 if(.not.rtdb_get(irtdb,'qmmm:nobq',mt_int,1,nobq)) 1155 + nobq=1 1156 uqmatm=cau2kj*uqmatm 1157 iquant=0 1158 else 1159 iquant=1 1160 endif 1161 endif 1162c 1163 if(.not.rtdb_cget(irtdb,'task:operation',1,operat)) 1164c + call md_abort('md_rtdbin: rtdg_cget failed',138) 1165 + operat='energy' 1166c 1167 if(operat.eq.'energy'.or.operat.eq.'ENERGY') then 1168 ntype=0 1169 endif 1170c 1171 if(operat.eq.'optimize'.or.operat.eq.'OPTIMIZE') then 1172 ntype=1 1173 if(nfrest.gt.0) nfqrs=nfrest 1174 nfvelo=0 1175 nfsvel=0 1176 endif 1177c 1178 if(operat.eq.'dynamics'.or.operat.eq.'DYNAMICS') then 1179 ntype=2 1180 mdtype=iset 1181 if(isetp1.eq.2) mdtype=4 1182 if(isetp1.eq.3) mdtype=5 1183 if(isetp1.eq.2.and.isetp2.eq.3) mdtype=6 1184 if(mdtype.gt.3) iset=1 1185 endif 1186c 1187 if(operat.eq.'thermodynamics'.or.operat.eq.'THERMODYNAMICS') then 1188 ntype=3 1189 endif 1190c 1191 ssshft=delta 1192c 1193 if(nodpme.gt.ngz) nodpme=ngz 1194 if(nodpme.gt.np) nodpme=np 1195 if(nodpme.eq.0) nodpme=np 1196 if(nodpme.gt.ngz) nodpme=ngz 1197c 1198c determine if this is a start, restart or continue 1199c 1200 call util_get_rtdb_state(irtdb,lstart,lrstrt,lcont) 1201c 1202 if(istart.gt.0) then 1203 lstart=istart.eq.0 1204 lrstrt=istart.eq.1 1205 lcont=istart.eq.2 1206 endif 1207c 1208 if(lstart) then 1209 nserie=0 1210 else 1211 if(lrstrt) then 1212 nserie=1 1213 else 1214 if(.not.lcont) 1215 + call md_abort('md_rtdbin: failed to determine rtdb state',0) 1216 nserie=2 1217 endif 1218 endif 1219c 1220c change lfnout if output to unit 6 is requested 1221c 1222 if(lfout6.ne.0) lfnout=6 1223c 1224c pairlist type 1225c 1226 lstype=1 1227 if(latom.eq.1) lstype=0 1228c 1229c Count radial distribution functions contributions 1230c ------------------------------------------------- 1231c 1232 ngc=0 1233 ngr=0 1234 ngrww=0 1235 ngrsw=0 1236 ngrss=0 1237 if(nfrdf.gt.0) then 1238 open(unit=lfnrdi,file=filrdi,form='formatted',status='old') 1239 rewind(unit=lfnrdi) 1240 1 continue 1241 read(lfnrdi,*,end=9999) igt,igr 1242 ngc=ngc+1 1243 if(igr.gt.ngr) ngr=igr 1244 if(igt.eq.1) ngrww=ngrww+1 1245 if(igt.eq.2) ngrsw=ngrsw+1 1246 if(igt.eq.3) ngrss=ngrss+1 1247 goto 1 1248 9999 continue 1249 close(unit=lfnrdi) 1250c 1251 if(ngc.gt.0) then 1252 if(ngl.eq.0) ngl=1000 1253 if(rrdf.lt.small) then 1254 if(ngrww.gt.0.and.rrdf.lt.rshort) rrdf=rshort 1255 if(ngrsw.gt.0.and.rrdf.lt.rshort) rrdf=rshort 1256 if(ngrss.gt.0.and.rrdf.lt.rshort) rrdf=rshort 1257 endif 1258 else 1259 ngl=1 1260 nfrdf=0 1261 endif 1262c 1263 endif 1264c routine md_chckin performs a check of selected input 1265c 1266c MCTI initializations 1267c 1268 if(mdtype.ge.7.and.ntype.eq.3.and.mrun.eq.0) mrun=maxlam 1269 if(mdtype.ge.7.and.ntype.eq.3.and.nffree.lt.0) nffree=1 1270 if(mdtype.ge.7.and.ntype.eq.3.and.nfoutp.lt.0) nfoutp=1 1271 if(mdtype.ge.7.and.ntype.eq.3.and.nfstat.lt.0) nfstat=1 1272 if(mdtype.ge.7.and.ntype.eq.3.and.ldacq.eq.0) ldacq=mdacq 1273 if(mdtype.ge.7.and.ntype.eq.3.and.macfl.le.mdacq) macfl=mdacq+1 1274 if(mdtype.ge.7.and.maxlam.le.1) 1275 + call md_abort('Number of MCTI ensembles too small:',maxlam) 1276c 1277c MD initializations 1278c 1279 if(nfoutp.lt.0) nfoutp=mdacq/100 1280 if(nfstat.lt.0) nfstat=mdacq/10 1281c 1282c square SHAKE tolerances 1283c 1284 tlwsha=tlwsha*tlwsha 1285 tlssha=tlssha*tlssha 1286c 1287c handle zero frequencies 1288c 1289 iint=0 1290 if(nfprop.lt.0) then 1291 iint=1 1292 nfprop=-nfprop 1293 endif 1294 if(nserie.ne.0) irr=1 1295c 1296 if(mrun.lt.0) then 1297 mrun=-mrun 1298 noone=mrun-1 1299 msplit=3 1300 endif 1301c 1302 if(nfdip.eq.0) then 1303 nfdip=1 1304 ndip=0 1305 endif 1306c 1307 if(nfkirk.eq.0) then 1308 nfkirk=1 1309 nkirk=0 1310 endif 1311c 1312 ivreas=nfgaus 1313 if(nfgaus.lt.0) nfgaus=-nfgaus 1314c 1315c only 3D-periodic boundary systems can be at constant pressure 1316c 1317c if(npbtyp.ne.1) ipscal=0 1318c 1319 lstate=rtdb_parallel(.true.) 1320c 1321c get variables from restart file if not initial run 1322c 1323 if(nserie.gt.0) call md_rdrest(lfnrst,filrst) 1324c 1325 endif 1326c 1327c broadcast to all nodes 1328c 1329 if(np.gt.1) then 1330 nbytes=loc(inp_ptr)-loc(ptol) 1331 call ga_brdcst(mrg_d36,ptol,nbytes,0) 1332 nbytes=7*ma_sizeof(mt_int,1,mt_byte) 1333 call ga_brdcst(mrg_d37,npx,nbytes,0) 1334 call util_char_ga_brdcst(mrg_d38,filnam,0) 1335 endif 1336c 1337 lesp=.false. 1338 lqmd=iquant.ne.0 1339 lqmmm=iqmmm.ne.0 1340 lpert2=ntype.eq.3.or.(iset.eq.1.and.(isetp1.eq.2.or.isetp2.eq.2)) 1341 lpert3=ntype.eq.3.or.(iset.eq.1.and.(isetp1.eq.3.or.isetp2.eq.3)) 1342c 1343 lpme=ipme.gt.0 1344 ltwin=rlong.gt.rshort.or.lpme 1345c 1346 lpack=nopack.eq.0 1347 if(lqmd) lpack=.false. 1348c 1349 msa=msareq 1350 mwm=mwmreq 1351c 1352 factgf=one 1353 if(iguide.gt.0) factgf=tstep/tguide 1354 factgg=one-factgf 1355c 1356 costio=zero 1357 lpmfc=.true. 1358c 1359 return 1360 end 1361 subroutine md_print() 1362c 1363 implicit none 1364c 1365#include "md_common.fh" 1366c 1367 character*22 ctype 1368c 1369 if(me.eq.0) then 1370 if(lfnout.ne.6) 1371 + open(unit=lfnout,file=filout(1:index(filout,' ')-1), 1372 + form='formatted',status='unknown') 1373c 1374 if(ntype.eq.0) then 1375 ctype='ENERGY EVALUATION ' 1376 elseif(ntype.eq.1) then 1377 ctype='ENERGY MINIMIZATION ' 1378 elseif(ntype.eq.2) then 1379 ctype='MOLECULAR DYNAMICS ' 1380 elseif(ntype.eq.3) then 1381 ctype='FREE ENERGY EVALUATION' 1382 else 1383 call md_abort('Unknown calculation type',ntype) 1384 endif 1385c 1386 call swatch(today,now) 1387 write(lfnout,1000) ctype,today,now 1388 1000 format(/,1x,a22,t110,2a10) 1389c 1390 write(lfnout,1001) titinp(1),datinp,timinp,titinp(2),titinp(3) 1391 1001 format(/,' Title ',t10,a,t110,2a10,/,t10,a,/,t10,a) 1392c 1393 write(lfnout,1002) filnam(1:index(filnam,' ')-1) 1394 1002 format(/,' System ',a) 1395c 1396 if(ntype.le.2) then 1397 write(lfnout,1003) iset 1398 1003 format(/,' Force field parameter set ',i4) 1399 endif 1400c 1401 if(nserie.eq.0) then 1402 write(lfnout,1026) 1403 1026 format(/,' Initial simulation',/) 1404 elseif(nserie.eq.1) then 1405 write(lfnout,1027) 1406 1027 format(/,' Restart simulation',/) 1407 else 1408 write(lfnout,1028) 1409 1028 format(/,' Continuation simulation',/) 1410 endif 1411c 1412 if(ntype.eq.3) then 1413 if(mropt.eq.0) then 1414 write(lfnout,1032) 1415 1032 format(' Multiple run initial simulation') 1416 elseif(mropt.eq.1) then 1417 write(lfnout,1033) 1418 1033 format(' Multiple run simulation with initial conditions') 1419 else 1420 write(lfnout,1034) 1421 1034 format(' Multiple run extension simulation') 1422 endif 1423 write(lfnout,1025) mrun 1424 1025 format(' Number of runs ',i7) 1425 endif 1426 if(ntype.ge.2) then 1427 write(lfnout,1024) mequi,mdacq 1428 1024 format(' Number of equilibration steps ',i7,/ 1429 + ' Number of data gathering steps',i7) 1430 endif 1431c 1432 if(ntype.eq.3.and.issscl.gt.0) then 1433 write(lfnout,1029) ssshft 1434 1029 format(/,' Separation shifted scaling, delta ',f12.6,' nm**2') 1435 endif 1436c 1437 if(ipme.gt.0) then 1438 write(lfnout,1030) morder,ngx,ngy,ngz,nodpme,pmetol 1439 1030 format(/,' Particle-mesh Ewald, spline to order ',i5,/, 1440 + ' Grid size ',i5,'x',i5,'x',i5,/, 1441 + ' Number of processors for p-FFT ',i5,/, 1442 + ' Tolerance at cutoff ',1pe12.5) 1443 if(imfft.le.1) then 1444 write(lfnout,1036) 1445 1036 format(' Using NWChem 3D-pFFT') 1446 else 1447 write(lfnout,1037) 1448#if defined(ESSL) 1449 1037 format(' Using PESSL 3D-pFFT') 1450#else 1451 1037 format(' Using system specific 3D-pFFT') 1452#endif 1453 endif 1454 endif 1455c 1456 if(ipolt.eq.1) then 1457 write(lfnout,1031) mpolit,ptol 1458 1031 format(/,' Polarization model, maximum iterations ',i5,/, 1459 + ' Tolerance ',1pe12.5) 1460 endif 1461c 1462 if(.not.ltwin) then 1463 write(lfnout,1004) rshort 1464 1004 format(/,' Cutoff radius ',f12.6, ' nm') 1465 else 1466 write(lfnout,1005) rshort,rlong 1467 1005 format(/,' Cutoff radius short range ',f12.6, ' nm', 1468 + /,' Cutoff radius long range ',f12.6, ' nm') 1469 endif 1470c 1471 if(ifield.gt.0) then 1472 write(lfnout,1035) xfield,xffreq,xfvect 1473 1035 format(/,' External electrostatic field strength ',f12.6,/, 1474 + 25x,'frequency ',f12.6,/, 1475 + 25x,'field vector ',3f12.6) 1476 endif 1477c 1478 if(ntype.ge.2) then 1479 if(ivreas.ne.0) then 1480 if(ivreas.gt.0) then 1481 write(lfnout,1022) tgauss,ivreas 1482 1022 format(/,' Velocity reassignment temperature ',f12.6,' K',/, 1483 + ' Velocity reassignment frequency ',i5,/) 1484 else 1485 write(lfnout,1122) tgauss 1486 1122 format(/,' Velocity reassignment temperature ',f12.6,' K',/, 1487 + ' Velocity reassignment in first step only',/) 1488 endif 1489 endif 1490 if(itscal.gt.0) then 1491 write(lfnout,1021) tmpext1,tmpext2,tmprlx,tmsrlx,tann1,tann2 1492 1021 format(/,' Isothermal ensemble external temperature ', 1493 + 2f12.6,' K',/, 1494 + ' Temperature relaxation time solvent ',f12.6,' ps',/, 1495 + ' Temperature relaxation time solute ',f12.6,' ps',/, 1496 + ' Temperature annealing between times ',2f12.6,' ps') 1497 endif 1498 if(ipscal.ne.0) then 1499 write(lfnout,1023) prsext,prsrlx,compr 1500 1023 format(' Isobaric ensemble external pressure ',1pe12.5,' Pa',/, 1501 + ' Pressure relaxation time ',0pf12.6,' ps',/, 1502 + ' Compressebility ',1pe12.5,' m**2/N') 1503 if(ipopt.eq.12) write(lfnout,2023) 1504 2023 format(' Pressure scaling in x and y only') 1505 if(ipopt.eq.3) write(lfnout,3023) 1506 3023 format(' Pressure scaling in z only') 1507 if(ipopt.eq.123) write(lfnout,3024) 1508 3024 format(' Pressure scaling in x and y separate from z') 1509 endif 1510 if(iguide.gt.0) then 1511 write(lfnout,1038) fguide,tguide 1512 1038 format(/,' Self-guided molecular dynamics with force coupling of', 1513 + f8.5,' and relaxation time of',f8.5,' ps') 1514 endif 1515 if(mdalgo.eq.1) then 1516 write(lfnout,1039) 1517 1039 format(/,' Leap-frog integration') 1518 else 1519 write(lfnout,1040) 1520 1040 format(/,' Leap-frog integration (Brown-Clarke)') 1521 endif 1522 endif 1523c 1524 if(scaleq.ge.zero) then 1525 write(lfnout,1041) scaleq 1526 1041 format(/,' Solute charges scaled by factor ',f8.5) 1527 endif 1528c 1529 if(ntype.eq.1) then 1530 if(msdit.gt.0) then 1531 write(lfnout,1006) msdit,dx0sd,dxmsd,dxsdmx 1532 1006 format(/,' Maximum number of steepest descent steps ',i5,/, 1533 + ' Initial step size ',f12.6,' nm',/, 1534 + ' Maximum step size ',f12.6,' nm',/, 1535 + ' Minimum step size ',f12.6,' nm') 1536 endif 1537 if(mcgit.gt.0) then 1538 write(lfnout,1007) mcgit,ncgcy,dx0cg,dxsdmx 1539 1007 format(/,' Maximum number of conjugate gradient steps ',i5,/, 1540 + ' Number of conjugate gradient steps per cycle ',i5,/, 1541 + ' Initial interval size ',f12.6,' nm',/, 1542 + ' Minimum step size ',f12.6,' nm') 1543 endif 1544 endif 1545c 1546 write(lfnout,1008) mshitw,tlwsha,mshits,tlssha 1547 1008 format(/,' Maximum number of solvent SHAKE iterations ',i5, 1548 + ', solvent SHAKE tolerance ',f12.6,' nm',/, 1549 + ' Maximum number of solute SHAKE iterations ',i5, 1550 + ', solute SHAKE tolerance ',f12.6,' nm',/) 1551c 1552 if(ntype.ge.2) then 1553 write(lfnout,1019) nfpair 1554 1019 format(' Frequency update pair lists ',t55,i5) 1555 if(ltwin) write(lfnout,1020) nflong 1556 1020 format(' Frequency update long range forces ',t55,i5) 1557 write(lfnout,1017) nfslow 1558 1017 format(' Frequency removal overall motion ',t55,i5) 1559 write(lfnout,1018) nfcntr 1560 1018 format(' Frequency solute centering ',t55,i5) 1561 if(icentr.eq.1) write(lfnout,1818) 1562 1818 format(' Centering in z-direction only') 1563 if(icentr.eq.2) write(lfnout,1819) 1564 1819 format(' Centering in xyplane only') 1565 write(lfnout,1009) nfoutp 1566 1009 format(' Frequency printing step information ',t55,i5) 1567 write(lfnout,1010) nfstat 1568 1010 format(' Frequency printing statistical information ',t55,i5) 1569 endif 1570 if(ntype.eq.1) then 1571 write(lfnout,1011) nfqrs 1572 1011 format(' Frequency recording minimum energy structure ',t55,i5) 1573 endif 1574 if(ntype.ge.2) then 1575 write(lfnout,1012) nfrest 1576 1012 format(' Frequency recording restart file ',t55,i5) 1577 write(lfnout,1013) nfcoor 1578 1013 format(' Frequency recording system coordinates ',t55,i5) 1579 write(lfnout,1014) nfscoo 1580 1014 format(' Frequency recording solute coordinates ',t55,i5) 1581 write(lfnout,1015) nfvelo 1582 1015 format(' Frequency recording system velocities ',t55,i5) 1583 write(lfnout,1016) nfsvel 1584 1016 format(' Frequency recording solute velocities ',t55,i5) 1585 write(lfnout,1115) nfforc 1586 1115 format(' Frequency recording system forces ',t55,i5) 1587 write(lfnout,1116) nfsfor 1588 1116 format(' Frequency recording solute forces ',t55,i5) 1589 endif 1590c 1591 write(lfnout,2000) 1592 2000 format(/,' LOAD BALANCING',/) 1593c 1594 1595 if(loadb.eq.0) write(lfnout,2001) 1596 2001 format(' None') 1597 if(loadb.eq.2) write(lfnout,2002) 1598 2002 format(' Redistribution of inter-processor box pairs') 1599 if(loadb.eq.1) write(lfnout,2003) factld 1600 2003 format(' Resizing of processor domains: smallest scaling ', 1601 + t55,f12.6) 1602 if(loadb.eq.3) write(lfnout,2004) lbpair,factld 1603 2004 format(' Redistribution of inter-processor box pairs: retry ', 1604 + t55,i5,/,' Resizing of processor domains: smallest scaling ', 1605 + t55,f12.6) 1606 if(isload.gt.0) write(lfnout,2104) 1607 2104 format(' Resizing in z-dimension only') 1608c 1609 if(ireset.gt.0) write(lfnout,2005) 1610 2005 format(/,' Load balancing reset') 1611c 1612 write(lfnout,2012) nfload 1613 2012 format(/,' Load balancing frequency ',i5) 1614c 1615 if(itload.eq.0) write(lfnout,2006) 1616 2006 format(/,' Load balancing based on last synchronization time') 1617 if(itload.eq.1) write(lfnout,2007) 1618 2007 format(/,' Load balancing based on minimum synchronization time') 1619 if(itload.eq.2) write(lfnout,2008) 1620 2008 format(/,' Load balancing based on average synchronization time') 1621 if(itload.eq.3) write(lfnout,2009) 1622 2009 format(/,' Load balancing based on average synchronization time', 1623 + ' with minimum synchronization time on processor 0') 1624c 1625 if(ioload.eq.1) write(lfnout,2010) 1626 2010 format(/,' Time for I/O counted at each step') 1627 if(ioload.eq.1) write(lfnout,2011) 1628 2011 format(/,' Experimental load balancing') 1629c 1630 if(nfhop.gt.0) then 1631c 1632 write(lfnout,3000) 1633 3000 format(/,' QHOP PROTON HOPPING',/) 1634c 1635 write(lfnout,3001) nfhop 1636 3001 format(' Frequency hopping attempts ',t55,i5) 1637 write(lfnout,3002) rhop 1638 3002 format(' Donor-acceptor cutoff distance',t55,f12.6,' nm') 1639 write(lfnout,3003) thop 1640 3003 format(' Minimum time before backhop',t55,f12.6,' ps') 1641c 1642 endif 1643c 1644 endif 1645c 1646 return 1647 end 1648 subroutine md_rdrest(lfn,fil) 1649c 1650 implicit none 1651c 1652#include "md_common.fh" 1653c 1654 integer lfn 1655 character*255 fil 1656c 1657 character*13 string 1658c 1659 if(me.eq.0) then 1660 open(unit=lfn,file=fil(1:index(fil,' ')-1), 1661 + status='old',form='formatted',err=9999) 1662 rewind(lfn) 1663c 1664 1 continue 1665 read(lfn,1000,end=9997) string 1666 1000 format(a13) 1667 if(string.ne.'restart input') goto 1 1668 read(lfn,1001,end=9998,err=9998) ntype,mdtype 1669 1001 format(11i7) 1670 read(lfn,1001) nfpair,nflong 1671 read(lfn,1001) lwtype,lstype,nfrest,keepr 1672 if(nserie.eq.1) then 1673 read(lfn,1002) krun,kequi,kdacq,mrun,mequi,mdacq,ldacq 1674 1002 format(7i7) 1675 else 1676 read(lfn,1002) krun,kequi,kdacq 1677 endif 1678 read(lfn,1003) stime,tstep 1679 1003 format(2f12.6) 1680 read(lfn,1004) rshort,rlong 1681 1004 format(2f12.6) 1682 read(lfn,1005) mshitw,tlwsha 1683 read(lfn,1005) mshits,tlssha 1684 1005 format(i7,f12.6) 1685 read(lfn,1006) ipscal,prsext,prsrlx,compr,ipopt 1686 1006 format(i5,e12.5,f12.6,e12.5,i5) 1687 read(lfn,1007) itscal,tmpext1,tmprlx,tmsrlx,tmpext2, 1688 + tann1,tann2 1689 1007 format(i5,6f12.6) 1690 read(lfn,1008) nfgaus,ivopt,tgauss,iseed 1691 1008 format(2i7,f12.6,i12) 1692 read(lfn,1009) nfoutp,nfstat,nfprop,nfnewf,ibatch 1693 1009 format(11i7) 1694 read(lfn,1009) ibinar,iformt 1695 read(lfn,1009) nfcoor,nfscoo,nfvelo,nfsvel,nfforc,nfsfor 1696 read(lfn,1009) nffree 1697 read(lfn,1009) nfcntr,nfslow 1698 read(lfn,1010) nfrdf,numrdf,ngc,ngr,ngl,ngrww,ngrsw,ngrss 1699 1010 format(8i7) 1700 read(lfn,1011) rrdf,drdf 1701 1011 format(2f12.6) 1702 read(lfn,1012) numdis,lendis 1703 1012 format(2i7) 1704 read(lfn,1013) numhis,lenhis 1705 1013 format(11i7) 1706 read(lfn,1014) nfdip,ndip,rdip 1707 1014 format(2i7,f12.6) 1708 read(lfn,1015) nfkirk,nkirk,rkirk 1709 1015 format(2i7,f12.6) 1710 endif 1711c 1712 9997 continue 1713 if(me.eq.0) close(unit=lfn) 1714 return 1715c 1716 9998 continue 1717 call md_abort('Unable to read restart file in md_rest ',me) 1718 return 1719 9999 continue 1720 call md_abort('Unable to open restart file in md_rest ',me) 1721 return 1722 end 1723 subroutine md_wtrest(lfn) 1724c 1725 implicit none 1726c 1727#include "md_common.fh" 1728c 1729 integer lfn 1730c 1731 if(me.ne.0) return 1732c 1733 write(lfn,1000) 1734 1000 format('restart input') 1735 write(lfn,1001) ntype,mdtype 1736 1001 format(11i7) 1737 write(lfn,1001) nfpair,nflong 1738 write(lfn,1001) lwtype,lstype,nfrest,keepr 1739 write(lfn,1002) irun,iequi,idacq,mrun,mequi,mdacq,ldacq 1740 1002 format(11i7) 1741 write(lfn,1003) stime,tstep 1742 1003 format(2f12.6) 1743 write(lfn,1004) rshort,rlong 1744 1004 format(2f12.6) 1745 write(lfn,1005) mshitw,tlwsha 1746 write(lfn,1005) mshits,tlssha 1747 1005 format(i7,f12.6) 1748 write(lfn,1006) ipscal,prsext,prsrlx,compr,ipopt 1749 1006 format(i5,e12.5,f12.6,e12.5,i5) 1750 write(lfn,1007) itscal,tmpext1,tmprlx,tmsrlx,tmpext2, 1751 + tann1,tann2 1752 1007 format(i5,6f12.6) 1753 tmpext=tmpext1 1754 write(lfn,1008) nfgaus,ivopt,tgauss,iseed 1755 1008 format(2i7,f12.6,i12) 1756 write(lfn,1009) nfoutp,nfstat,nfprop,nfnewf,ibatch 1757 1009 format(11i7) 1758 write(lfn,1009) ibinar,iformt 1759 write(lfn,1009) nfcoor,nfscoo,nfvelo,nfsvel,nfforc,nfsfor 1760 write(lfn,1009) nffree 1761 write(lfn,1009) nfcntr,nfslow 1762 write(lfn,1010) nfrdf,numrdf,ngc,ngr,ngl,ngrww,ngrsw,ngrss 1763 1010 format(8i7) 1764 write(lfn,1011) rrdf,drdf 1765 1011 format(2f12.6) 1766 write(lfn,1012) numdis,lendis 1767 1012 format(2i7) 1768 write(lfn,1013) numhis,lenhis 1769 1013 format(11i7) 1770 write(lfn,1014) nfdip,ndip,rdip 1771 1014 format(2i7,f12.6) 1772 write(lfn,1015) nfkirk,nkirk,rkirk 1773 1015 format(2i7,f12.6) 1774c 1775 return 1776 end 1777 subroutine md_fopen(lclose) 1778c 1779 implicit none 1780c 1781#include "md_common.fh" 1782c 1783 logical lclose 1784c 1785 if(nfnewf.gt.0) then 1786 ibatch=ibatch+1 1787 write(filtrj,1000) root(1:index(root,' ')-1), 1788 + ibatch,'.trj' 1789 write(filhop,1000) root(1:index(root,' ')-1), 1790 + ibatch,'.hop' 1791 write(filprp,1000) root(1:index(root,' ')-1), 1792 + ibatch,'.prp' 1793 write(filpmf,1000) root(1:index(root,' ')-1), 1794 + ibatch,'.pmf' 1795 write(filtim,1000) root(1:index(root,' ')-1), 1796 + ibatch,'.tim' 1797 write(rfile,1000) root(1:index(root,' ')-1), 1798 + ibatch,'.rst' 1799 1000 format(a,i3.3,a) 1800 endif 1801c 1802 if(ntype.ne.3) then 1803 if(nfcoor.gt.0.or.nfscoo.gt.0.or.nfvelo.gt.0.or.nfsvel.gt.0) then 1804 if(lclose) close(lfntrj) 1805 open(unit=lfntrj,file=filtrj(1:index(filtrj,' ')-1), 1806 + form='formatted',status='unknown') 1807 call cf_trjhdr(lfntrj) 1808 endif 1809 if(nfhop.gt.0) then 1810 if(lclose) close(lfnhop) 1811 open(unit=lfnhop,file=filhop(1:index(filhop,' ')-1), 1812 + form='formatted',status='unknown') 1813 endif 1814 if(nfprop.gt.0) then 1815 if(lclose) close(lfnprp) 1816 open(unit=lfnprp,file=filprp(1:index(filprp,' ')-1), 1817 + form='formatted',status='unknown') 1818 call prp_header() 1819 endif 1820 endif 1821 if(ntype.eq.2.and.iprpmf.ne.0) then 1822 if(lclose) close(lfnpmf) 1823 open(unit=iabs(lfnpmf),file=filpmf(1:index(filpmf,' ')-1), 1824 + form='formatted',status='unknown') 1825 endif 1826 if(nftime.gt.0) then 1827 if(lclose) close(lfntim) 1828 open(unit=lfntim,file=filtim(1:index(filtim,' ')-1), 1829 + form='formatted',status='unknown') 1830 call md_hdrtim() 1831 endif 1832c 1833 return 1834 end 1835 subroutine md_hdrtim() 1836c 1837 implicit none 1838c 1839#include "md_common.fh" 1840c 1841 write(lfntim,1200) np,56 1842 1200 format(2i5) 1843 write(lfntim,1201) 1844 write(lfntim,1202) 1845 write(lfntim,1203) 1846 write(lfntim,1204) 1847 write(lfntim,1205) 1848 write(lfntim,1206) 1849 write(lfntim,1207) 1850 write(lfntim,1208) 1851 write(lfntim,1209) 1852 write(lfntim,1210) 1853 write(lfntim,1211) 1854 write(lfntim,1212) 1855 write(lfntim,1213) 1856 write(lfntim,1214) 1857 write(lfntim,1215) 1858 write(lfntim,1216) 1859 write(lfntim,1217) 1860 write(lfntim,1218) 1861 write(lfntim,1219) 1862 write(lfntim,1220) 1863 write(lfntim,1221) 1864 write(lfntim,1222) 1865 write(lfntim,1223) 1866 write(lfntim,1224) 1867 write(lfntim,1225) 1868 write(lfntim,1226) 1869 write(lfntim,1227) 1870 write(lfntim,1228) 1871 write(lfntim,1229) 1872 write(lfntim,1230) 1873 write(lfntim,1231) 1874 write(lfntim,1232) 1875 write(lfntim,1233) 1876 write(lfntim,1234) 1877 write(lfntim,1235) 1878 write(lfntim,1236) 1879 write(lfntim,1237) 1880 write(lfntim,1238) 1881 write(lfntim,1239) 1882 write(lfntim,1240) 1883 write(lfntim,1241) 1884 write(lfntim,1242) 1885 write(lfntim,1243) 1886 write(lfntim,1244) 1887 write(lfntim,1245) 1888 write(lfntim,1246) 1889 write(lfntim,1247) 1890 write(lfntim,1248) 1891 write(lfntim,1249) 1892 write(lfntim,1250) 1893 write(lfntim,1251) 1894 write(lfntim,1252) 1895 write(lfntim,1253) 1896 write(lfntim,1254) 1897 write(lfntim,1255) 1898 write(lfntim,1256) 1899 1201 format(' 1 : Initialization velocities') 1900 1202 format(' 2 : Periodic Boundary Conditions') 1901 1203 format(' 3 : Dynamic Load Balancing') 1902 1204 format(' 4 : Atom Redistribution') 1903 1205 format(' 5 : Center of Mass Coordinates') 1904 1206 format(' 6 : Recording Trajectory') 1905 1207 format(' 7 : CAFE Initialization') 1906 1208 format(' 8 : QMD Forces Evaluation') 1907 1209 format(' 9 : Synchronization') 1908 1210 format(' 10 : SPACE Initialization') 1909 1211 format(' 11 : Synchronization') 1910 1212 format(' 12 : Collapse Option') 1911 1213 format(' 13 : External Fields') 1912 1214 format(' 14 : Multiprocs Data Initialization') 1913 1215 format(' 15 : Multiprocs Data Communication') 1914 1216 format(' 16 : Multiprocs Data Forces Evaluation') 1915 1217 format(' 17 : Restraints Data Initialization') 1916 1218 format(' 18 : Restraints Data Communication') 1917 1219 format(' 19 : Restraints Data Forces Evaluation') 1918 1220 format(' 20 : PMF Data Initialization') 1919 1221 format(' 21 : PMF Data Communication') 1920 1222 format(' 22 : PMF Data Forces Evaluation') 1921 1223 format(' 23 : Induced Dipoles Evaluation') 1922 1224 format(' 24 : PME Initialization') 1923 1225 format(' 25 : PME Charge Grid Evaluation') 1924 1226 format(' 26 : PME Charge Grid Communication') 1925 1227 format(' 27 : PME Charge Grid Retrieval') 1926 1228 format(' 28 : PME Reverse 3D-pFFT') 1927 1229 format(' 29 : PME Energy Evaluation') 1928 1230 format(' 30 : PME Synchronization') 1929 1231 format(' 31 : PME Forward 3D-pFFT') 1930 1232 format(' 32 : PME Synchronization') 1931 1233 format(' 33 : Local Coordinates Retrieval') 1932 1234 format(' 34 : Local Forces Evaluation') 1933 1235 format(' 35 : Local Force Accumulation') 1934 1236 format(' 36 : Non-Local Coordinates Retrieval') 1935 1237 format(' 37 : Non-Local Forces Evaluation') 1936 1238 format(' 38 : Non-Local Force Accumulation') 1937 1239 format(' 39 : PME Wait') 1938 1240 format(' 40 : PME Forces Evaluation') 1939 1241 format(' 41 : PME Forces Communication') 1940 1242 format(' 42 : PME Node Synchronization') 1941 1243 format(' 43 : PME Flag') 1942 1244 format(' 44 : Synchronization') 1943 1245 format(' 45 : SPACE Finalization') 1944 1246 format(' 46 : QM/MM Forces Evaluation') 1945 1247 format(' 47 : Forces Normalization') 1946 1248 format(' 48 : Guided MD Forces') 1947 1249 format(' 49 : MD Time Step') 1948 1250 format(' 50 : SHAKE') 1949 1251 format(' 51 : CAFE Finalization') 1950 1252 format(' 52 : Center of Mass Motion Removal') 1951 1253 format(' 53 : SPACE Update') 1952 1254 format(' 54 : Property Evaluation') 1953 1255 format(' 55 : Restart File Update') 1954 1256 format(' 56 : TOTAL TIME PER STEP') 1955c 1956 return 1957 end 1958 subroutine md_partition 1959c 1960 implicit none 1961c 1962#include "md_common.fh" 1963#include "global.fh" 1964#include "util.fh" 1965c size should equal maximum number of processors that could be in 1966c a group 1967 integer MAX_GRP_SIZE 1968 parameter (MAX_GRP_SIZE = 1000) 1969 integer grp_size, nproc 1970 integer list(MAX_GRP_SIZE) 1971 integer i, grp_handle 1972 real*8 rantmp 1973c 1974c write(*,1000) npg 1975c 1000 format(' The number of processor groups requested is ',i5) 1976c write(*,1001) ga_pgroup_get_world() 1977c 1001 format(' The world group handle is ',i5) 1978c write(*,1002) ga_pgroup_get_default() 1979c 1002 format(' The default group handle is ',i5) 1980c 1981 nproc = ga_nnodes() 1982 me = ga_nodeid() 1983 rantmp=util_random(iseed*(me+1)) 1984 if (mod(nproc,npg).ne.0) then 1985 call md_abort('Number of processor groups is not a multiple', 1986 + ' of number of processors',ga_nodeid()) 1987 endif 1988 grp_size = nproc/npg 1989 ipg = (me - mod(me,grp_size))/grp_size 1990 do i = 1, grp_size 1991 list(i) =ipg*grp_size + i - 1 1992 end do 1993 grp_handle = ga_pgroup_create(list,grp_size) 1994 call ga_pgroup_set_default(grp_handle) 1995 np = ga_nnodes() 1996 me = ga_nodeid() 1997 return 1998 end 1999c 2000 subroutine md_partition_end 2001c 2002 implicit none 2003c 2004#include "md_common.fh" 2005#include "global.fh" 2006 integer grp_handle 2007 logical status 2008c 2009 grp_handle = ga_pgroup_get_default() 2010 call ga_pgroup_set_default(ga_pgroup_get_world()) 2011 status = ga_pgroup_destroy(grp_handle) 2012c 2013 call ga_sync() 2014c 2015 np = ga_nnodes() 2016 me = ga_nodeid() 2017 return 2018 end 2019