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