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