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