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