1c
2c $Id$
3c
4
5      SUBROUTINE md_inpt(inputfile,
6     $                   iseed,tstep,nstep,nequl,nprnt,ntype,ncons,
7     $    consatm,nbond,bondatm,nshel,shelatm,lveloc)
8
9      implicit none
10      character*(*) inputfile
11
12      include 'p_input.inc'
13      include 'p_array.inc'
14      include 'cm_atom.inc'
15      include 'cm_latt.inc'
16      include 'cm_temp.inc'
17      include 'cm_ewld.inc'
18      include 'cm_cuto.inc'
19      include 'cm_pote.inc'
20      include 'cm_cons.inc'
21      include 'cm_bond.inc'
22      include 'cm_shel.inc'
23
24      integer i,j,k
25      integer ierr,ifield,istart,inum_char,infile
26      integer iseed,ntype,nequl,nstep,nprnt
27      integer ncons,consatm,nbond,bondatm,nshel,shelatm
28      integer npote,potindex,pkey,patm1,patm2,imin,imax,idif
29
30      real*8 tstep,para1,para2,para3
31
32      character word*4,irecord*1,filename*20
33
34      logical lend,lcoord,lveloc,lvectr
35
36      dimension consatm(mxcons,2),bondatm(mxbond,2),shelatm(mxshel,2)
37      dimension istart(irec_len),inum_char(irec_len),irecord(irec_len)
38
39      open(unit=input,file=inputfile)
40      infile=input
41
42      nstep=0
43      nequl=0
44      nprnt=1
45      ntype=0
46      npote=0
47      ncons=0
48      nbond=0
49      nshel=0
50      natms=0
51      iseed=620419483
52      kmaxx=0
53      kmaxy=0
54      kmaxz=0
55      tstep=0.0
56      temp=0.0
57      rcut=0.0
58      vcut=0.0
59      alpha=0.0
60      lcoord=.false.
61      lveloc=.false.
62      lvectr=.false.
63
64!##########################################################################
65 100  call read_lin(infile,ifield,istart,inum_char,irecord,lend)
66      if(lend)goto 200
67      if(ifield.gt.0)then
68       call ex_4char(istart(1),irecord,word)
69!##########################################################################
70       if(word.eq.'INCL')then
71        j=0
72        do k=istart(2),istart(2)+inum_char(2)
73         j=j+1
74         filename(j:j)=irecord(k)
75        enddo
76        write(output,"(/,1x,'Openning include file :',a20)")
77     $   filename(1:j)
78        open(unit=iinclude,file=filename(1:j))
79        infile=iinclude
80        goto 100
81!##########################################################################
82       elseif(word.eq.'VECT')then
83        lvectr=.true.
84        read(infile,*)latt(1,1),latt(2,1),latt(3,1)
85        read(infile,*)latt(1,2),latt(2,2),latt(3,2)
86        read(infile,*)latt(1,3),latt(2,3),latt(3,3)
87        write(output,"(/,1x,'Simulation cell vectors')")
88        write(output,'(3g20.10)')latt(1,1),latt(2,1),latt(3,1)
89        write(output,'(3g20.10)')latt(1,2),latt(2,2),latt(3,2)
90        write(output,'(3g20.10)')latt(1,3),latt(2,3),latt(3,3)
91!##########################################################################
92       elseif(word.eq.'COOR')then
93        call ex_inter(inum_char(2),istart(2),irecord,natms,ierr)
94        lcoord=.true.
95        if(ierr.eq.1)then
96         write(output,*)'Invalid number of particles'
97        endif
98        if(natms.gt.mxatms)then
99         write(output,*)'Increase mxatms to ',natms
100         stop
101        endif
102        do j=1,natms
103         read(infile,*,end=300)atmtype(j),ccc(j,1),ccc(j,2),ccc(j,3)
104        enddo
105        write(output,"(/,1x,i9,' atom coordinates found')")natms
106!##########################################################################
107       elseif(word.eq.'VELO')then
108        call ex_inter(inum_char(2),istart(2),irecord,natms,ierr)
109        lveloc=.true.
110        if(ierr.eq.1)then
111         write(output,*)'Invalid number of particles'
112         stop
113        endif
114        if(natms.gt.mxatms)then
115         write(output,*)'Increase mxatms to ',natms
116         stop
117        endif
118        do j=1,natms
119         read(infile,*,end=300)vvv(j,1),vvv(j,2),vvv(j,3)
120        enddo
121        write(output,"(/,1x,i9,' atom velocities found')")natms
122!##########################################################################
123       elseif(word.eq.'ATOM')then
124        call ex_inter(inum_char(2),istart(2),irecord,ntype,ierr)
125        if(ierr.eq.1)then
126         write(output,*)'Invalid number of particle types'
127         stop
128        endif
129        if(ntype.gt.mxtype)then
130         write(output,*)'Increase mxtype to ',ntype
131         stop
132        endif
133        write(output,"(/,1x,'Atom types: ',i6)")ntype
134        potindex=0
135        do j=1,ntype
136         read(infile,*,end=300)typname(j),typmass(j),typchge(j),
137     $    typmol(j)
138         write(output,'(a8,1x,2(g20.10,1x),i6)')typname(j),typmass(j),
139     $      typchge(j),typmol(j)
140         do k=j,ntype
141          potindex=potindex+1
142          potkey(potindex)=0
143         enddo
144        enddo
145!##########################################################################
146       elseif(word.eq.'POTE')then
147        call ex_inter(inum_char(2),istart(2),irecord,npote,ierr)
148        if(npote.gt.mxpote)then
149         write(output,*)'Increase mxpote to ',npote
150         stop
151        endif
152        if(ntype.eq.0)then
153         write(output,*)'Atom type info needs to be provided first'
154         stop
155        endif
156        write(output,"(/,1x,'Potential parameters:')")
157        do j=1,npote
158         read(infile,*,end=300)pkey,patm1,patm2,para1,para2,para3
159         imin=min(patm1,patm2)
160         imax=max(patm1,patm2)
161         idif=imax-imin
162         potindex=0
163         do k=1,imin-1
164          potindex=potindex+(ntype-k+1)
165         enddo
166         potindex=potindex+idif+1
167         potkey(potindex)=pkey
168         potpar(potindex,1)=para1
169         potpar(potindex,2)=para2
170         potpar(potindex,3)=para3
171        enddo
172        potindex=0
173        do j=1,ntype
174         do k=j,ntype
175         potindex=potindex+1
176         if(potkey(potindex).gt.0)then
177          write(output,'(3(i3,1x),3(1x,g20.10))')j,k,potkey(potindex),
178     $     potpar(potindex,1),potpar(potindex,2),potpar(potindex,3)
179          endif
180         enddo
181        enddo
182!##########################################################################
183       elseif(word.eq.'CONS')then
184        call ex_inter(inum_char(2),istart(2),irecord,ncons,ierr)
185        if(ierr.eq.1)then
186         write(output,*)'Invalid number of constraints'
187         stop
188        endif
189        if(ncons.gt.mxcons)then
190         write(output,*)'Increase mxcons to ',ncons
191         stop
192        endif
193        write(output,"(/,1x,'Distance constraints: ',i6)")ncons
194        do j=1,ncons
195         read(infile,*,end=300)consatm(j,1),consatm(j,2),consdist(j)
196         write(output,'(2(i6,1x),g20.10)')consatm(j,1),consatm(j,2),
197     $      consdist(j)
198        enddo
199!##########################################################################
200       elseif(word.eq.'BOND')then
201        call ex_inter(inum_char(2),istart(2),irecord,nbond,ierr)
202        if(ierr.eq.1)then
203         write(output,*)'Invalid number of bonds'
204         stop
205        endif
206        if(nbond.gt.mxbond)then
207         write(output,*)'Increase mxbond to ',nbond
208         stop
209        endif
210        write(output,"(/,1x,'Bonded interactions: ',i6)")nbond
211        do j=1,nbond
212         read(infile,*,end=300)bondatm(j,1),bondatm(j,2),bondfrce(j),
213     $      bonddist(j)
214         write(output,'(2(i6,1x),2(g20.10,1x))')bondatm(j,1),
215     $      bondatm(j,2),bondfrce(j),bonddist(j)
216        enddo
217!##########################################################################
218       elseif(word.eq.'SHEL')then
219        call ex_inter(inum_char(2),istart(2),irecord,nshel,ierr)
220        if(ierr.eq.1)then
221         write(output,*)'Invalid number of shells'
222         stop
223        endif
224        if(nshel.gt.mxshel)then
225         write(output,*)'Increase mxshel to ',nshel
226         stop
227        endif
228        write(output,"(/,1x,'Number of core-shell units: ',i6)")nshel
229        do j=1,nshel
230         read(infile,*,end=300)shelatm(j,1),shelatm(j,2),shelfrce(j)
231         write(output,'(2(i6,1x),g20.10)')shelatm(j,1),shelatm(j,2),
232     $      shelfrce(j)
233        enddo
234!##########################################################################
235       elseif(word.eq.'SEED')then
236        call ex_inter(inum_char(2),istart(2),irecord,iseed,ierr)
237        if(ierr.eq.1)then
238         write(output,*)'Invalid seed for random number generator'
239         stop
240        endif
241        write(output,"(/,1x,'Seed for random number generator :',i18)")
242     $   iseed
243!##########################################################################
244       elseif(word.eq.'TEMP')then
245        call ex_1real(inum_char(2),istart(2),irecord,temp,ierr)
246        if(ierr.eq.2)then
247         write(output,*)'Invalid temperature'
248         stop
249        endif
250        write(output,"(/,1x,'Temperature :',3x,f12.4)")temp
251!##########################################################################
252       elseif(word.eq.'CUTO')then
253        call ex_1real(inum_char(2),istart(2),irecord,rcut,ierr)
254        if(ierr.eq.2)then
255         write(output,*)'Invalid cutoff distance'
256         stop
257        endif
258        write(output,"(/,1x,'Cutoff distance :',3x,f12.4)")rcut
259        rcutsq=rcut**2
260!##########################################################################
261       elseif(word.eq.'VERL')then
262        call ex_1real(inum_char(2),istart(2),irecord,vcut,ierr)
263        if(ierr.eq.2)then
264         write(output,*)'Invalid verlet shell cutoff'
265         stop
266        endif
267        write(output,"(/,1x,'Verlet shell cutoff :',3x,f12.4)")vcut
268        vcutsq=vcut**2
269!##########################################################################
270       elseif(word.eq.'EWAL')then
271        call ex_1real(inum_char(2),istart(2),irecord,alpha,ierr)
272        if(ierr.eq.2)then
273         write(output,*)'Invalid ewald parameter'
274         stop
275        endif
276        write(output,"(/,1x,'Ewald parameter :',3x,f12.4)")alpha
277!##########################################################################
278       elseif(word.eq.'KVEC')then
279        call ex_inter(inum_char(2),istart(2),irecord,kmaxx,ierr)
280        call ex_inter(inum_char(3),istart(3),irecord,kmaxy,ierr)
281        call ex_inter(inum_char(4),istart(4),irecord,kmaxz,ierr)
282        if(ierr.eq.1)then
283         write(output,*)'Invalid number of k-vectors'
284         stop
285        endif
286        if(kmaxx.gt.mxkvect)then
287         write(output,*)'Increase mxkvect to ',kmaxx
288         stop
289        endif
290        if(kmaxy.gt.mxkvect)then
291         write(output,*)'Increase mxkvect to ',kmaxy
292         stop
293        endif
294        if(kmaxz.gt.mxkvect)then
295         write(output,*)'Increase mxkvect to ',kmaxz
296         stop
297        endif
298        write(output,"(/,1x,'k-vectors :',3i4)")kmaxx,kmaxy,kmaxz
299!##########################################################################
300       elseif(word.eq.'STEP')then
301        call ex_inter(inum_char(2),istart(2),irecord,nstep,ierr)
302        if(ierr.eq.1)then
303         write(output,*)'Invalid number of MD steps'
304         stop
305        endif
306        write(output,"(/,1x,'Number of MD steps :',i12)")nstep
307!##########################################################################
308       elseif(word.eq.'EQUI')then
309        call ex_inter(inum_char(2),istart(2),irecord,nequl,ierr)
310        if(ierr.eq.1)then
311         write(output,*)'Invalid number of equilibration steps'
312         stop
313        endif
314        write(output,"(/,1x,'Number of equilibration steps :',i12)")
315     $ nequl
316!##########################################################################
317       elseif(word.eq.'PRIN')then
318        call ex_inter(inum_char(2),istart(2),irecord,nprnt,ierr)
319        if(ierr.eq.1)then
320         write(output,*)'Invalid printing frequency'
321         stop
322        endif
323        write(output,"(/,1x,'Print every ',i9,' steps')")nprnt
324!##########################################################################
325       elseif(word.eq.'TIME')then
326        call ex_1real(inum_char(2),istart(2),irecord,tstep,ierr)
327        if(ierr.eq.2)then
328         write(output,*)'Invalid integration timestep'
329         stop
330        endif
331        write(output,"(/,1x,'Integration timestep :',3x,f12.4)")tstep
332!##########################################################################
333       endif
334      endif
335
336      goto 100
337
338 200  continue
339
340      if(infile.eq.input)then
341       if(nstep.eq.0)then
342        write(output,*)'Number of MD steps has not been defined'
343        stop
344       elseif(ntype.eq.0)then
345        write(output,*)'Atom types have not been defined'
346        stop
347       elseif(natms.eq.0)then
348        write(output,*)'Number of atoms has not been defined'
349        stop
350       elseif(kmaxx.eq.0)then
351        write(output,*)'Number of k-vectors in x has not been defined'
352        stop
353       elseif(kmaxy.eq.0)then
354        write(output,*)'Number of k-vectors in y has not been defined'
355        stop
356       elseif(kmaxz.eq.0)then
357        write(output,*)'Number of k-vectors in z has not been defined'
358        stop
359       elseif(tstep.eq.0.0)then
360        write(output,*)'Integration timestep has not been defined'
361        stop
362       elseif(temp.eq.0.0)then
363        write(output,*)'Temperature has not been defined'
364        stop
365       elseif(rcut.eq.0.0)then
366        write(output,*)'Short-range cutoff has not been defined'
367        stop
368       elseif(vcut.eq.0.0)then
369        write(output,*)'Verlet list cutoff has not been defined'
370        stop
371       elseif(alpha.eq.0.0)then
372        write(output,*)'Ewald parameter has not been defined'
373        stop
374       elseif(.not.lcoord)then
375        write(output,*)'Could not find atomic coordinates'
376        stop
377       elseif(.not.lvectr)then
378        write(output,*)'Could not find simulation cell vectors'
379        stop
380       endif
381       write(output,*)
382       return
383      elseif(infile.eq.iinclude)then
384       infile=input
385       goto 100
386      endif
387
388 300  continue
389      write(output,"(/,1x,'EOF reached in keyword ',a4)")word
390
391      END
392