1c X. W. Zhou, xzhou@sandia.gov 2 include 'createAtoms.h' 3 call latgen() 4 call delete() 5 call hitvalue() 6 call disturb() 7 call systemshift() 8 call colect() 9 if (ilatseed .ne. 0) call randomize() 10 call velgen() 11 call outputfiles() 12 if (float(nntype(1)) .ne. 0.0) 13 * write(6,*)(float(nntype(i))/float(nntype(1)),i=1,ntypes) 14 write(6,*)'atoms of different species ',(nntype(i),i=1,ntypes) 15 stop 16 end 17cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 18 subroutine latgen() 19 include 'createAtoms.h' 20 character*8 lattype 21 common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3), 22 * periodicity(3),xbound(2),ybound(2),zbound(2),strain(3), 23 * delx(3),rcell(3),ccell(nelmax) 24 common /iran/ iseed 25 namelist /maincard/ ntypes,amass,ielement,ilatseed, 26 * perlb,perub,iseed,xy,xz,yz 27 namelist /latcard/ lattype,alat,xrot,yrot,zrot, 28 * periodicity,xbound,ybound,zbound, 29 * strain,delx 30 iseed=3245566 31 natoms=0 32 ilatseed=0 33 xy=1.0e12 34 xz=1.0e12 35 yz=1.0e12 36 read(5,maincard) 37 do i = 1,3 38 perlen(i)=perub(i)-perlb(i) 39 enddo 40 do i=1,ntypes 41 nntype(i)=0 42 enddo 431 continue 44 strain(1)=0.0 45 strain(2)=0.0 46 strain(3)=0.0 47 xbound(1)=0.0 48 xbound(2)=0.0 49 ybound(1)=0.0 50 ybound(2)=0.0 51 zbound(1)=0.0 52 zbound(2)=0.0 53 delx(1)=0.0 54 delx(2)=0.0 55 delx(3)=0.0 56 lattype='none' 57 read(5,latcard) 58 if (lattype .eq. 'none') goto 2 59 call latgen1() 60 call defvalue() 61 goto 1 622 continue 63 return 64 end 65cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 66 subroutine latgen1() 67 include 'createAtoms.h' 68 character*8 lattype 69 common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3), 70 * periodicity(3),xbound(2),ybound(2),zbound(2),strain(3), 71 * delx(3),rcell(3),ccell(nelmax) 72 common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3), 73 * avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100) 74 data xold/1.0,0.0,0.0/,yold/0.0,1.0,0.0/,zold/0.0,0.0,1.0/ 75 namelist /subcard/ rcell,ccell 76 77 call storelat(lattype,avec,bvec,cvec) 78 if (xbound(1) .eq. xbound(2)) then 79 xbound(1)=perlb(1) 80 xbound(2)=perub(1) 81 endif 82 if (ybound(1) .eq. ybound(2)) then 83 ybound(1)=perlb(2) 84 ybound(2)=perub(2) 85 endif 86 if (zbound(1) .eq. zbound(2)) then 87 zbound(1)=perlb(3) 88 zbound(2)=perub(3) 89 endif 90 vol=alat(1)*alat(2)*alat(3) 91 xnorm=sqrt((xrot(1)*alat(2)*alat(3))**2+ 92 * (xrot(2)*alat(1)*alat(3))**2+ 93 * (xrot(3)*alat(1)*alat(2))**2) 94 ynorm=sqrt((yrot(1)*alat(2)*alat(3))**2+ 95 * (yrot(2)*alat(1)*alat(3))**2+ 96 * (yrot(3)*alat(1)*alat(2))**2) 97 znorm=sqrt((zrot(1)*alat(2)*alat(3))**2+ 98 * (zrot(2)*alat(1)*alat(3))**2+ 99 * (zrot(3)*alat(1)*alat(2))**2) 100 periodicity(1)=vol*periodicity(1)/xnorm 101 periodicity(2)=vol*periodicity(2)/ynorm 102 periodicity(3)=vol*periodicity(3)/znorm 103 do 5 i=1,3 104 xrot(i)=vol*xrot(i)/(alat(i)*xnorm) 105 yrot(i)=vol*yrot(i)/(alat(i)*ynorm) 106 zrot(i)=vol*zrot(i)/(alat(i)*znorm) 1075 continue 108 do 10 i=1,3 109 do 10 j=1,3 110 roter(i,j)=xold(i)*xrot(j)+yold(i)*yrot(j)+zold(i)*zrot(j) 11110 continue 112 do 20 i=1,3 113 avecp(i)=0.0 114 bvecp(i)=0.0 115 cvecp(i)=0.0 116 do 20 j=1,3 117 avecp(i)=avecp(i)+roter(i,j)*avec(j)*0.5*alat(j) 118 bvecp(i)=bvecp(i)+roter(i,j)*bvec(j)*0.5*alat(j) 119 cvecp(i)=cvecp(i)+roter(i,j)*cvec(j)*0.5*alat(j) 12020 continue 121 do 30 i=1,3 122 avec(i)=avecp(i) 123 bvec(i)=bvecp(i) 124 cvec(i)=cvecp(i) 12530 continue 126 call applystrain() 127 nrange=20 128 nsmall=0 129 do 300 ii=-nrange,nrange 130 do 300 jj=-nrange,nrange 131 do 300 kk=-nrange,nrange 132 xt=ii*avec(1)+jj*bvec(1)+kk*cvec(1) 133 if (xt .ge. periodicity(1)-0.1 .or. xt .lt. -0.1) goto 300 134 yt=ii*avec(2)+jj*bvec(2)+kk*cvec(2) 135 if (yt .ge. periodicity(2)-0.1 .or. yt .lt. -0.1) goto 300 136 zt=ii*avec(3)+jj*bvec(3)+kk*cvec(3) 137 if (zt .ge. periodicity(3)-0.1 .or. zt .lt. -0.1) goto 300 138 nsmall=nsmall+1 139 rvs(1,nsmall)=xt 140 rvs(2,nsmall)=yt 141 rvs(3,nsmall)=zt 142300 continue 143 nxa=xbound(1)/periodicity(1)-2 144 nxb=xbound(2)/periodicity(1)+2 145 nya=ybound(1)/periodicity(2)-2 146 nyb=ybound(2)/periodicity(2)+2 147 nza=zbound(1)/periodicity(3)-2 148 nzb=zbound(2)/periodicity(3)+2 1491 rcell(1)=-100.0 150 read(5,subcard) 151 if (rcell(1) .eq. -100.0) then 152 return 153 endif 154 do 190 m = 2,ntypes 155 ccell(m)=min(1.0,ccell(m-1)+ccell(m)) 156190 continue 157 do 800 i=1,3 158 rcellp(i)=0.0 159 do 800 j=1,3 160 rcellp(i)=rcellp(i)+roter(i,j)*rcell(j)*alat(j) 161800 continue 162 do 801 i=1,3 163 rcell(i)=rcellp(i)*(1.0+strain(i)) 164801 continue 165 do 250 nn=1,nsmall 166 do 251 ii=nxa,nxb 167 xt=ii*periodicity(1)+rvs(1,nn)+rcell(1) 168 if (xt .ge. xbound(2) .or. xt .lt. xbound(1)) goto 251 169 do 252 jj=nya,nyb 170 yt=jj*periodicity(2)+rvs(2,nn)+rcell(2) 171 if (yt .ge. ybound(2) .or. yt .lt. ybound(1)) goto 252 172 do 253 kk=nza,nzb 173 zt=kk*periodicity(3)+rvs(3,nn)+rcell(3) 174 if (zt .ge. zbound(2) .or. zt .lt. zbound(1)) goto 253 175 tst=ranl() 176 it=0 177 do 201 m=ntypes,1,-1 178 if (tst .le. ccell(m)) it=m 179201 continue 180 natoms=natoms+1 181 if (natoms .gt. natmax) then 182 write(6,9601)natmax 1839601 format('number of atoms exceeds maximum dimension: ',i10) 184 endif 185 itype(natoms)=it 186 nhittag(natoms)=0 187 nntype(it)=nntype(it)+1 188 rv(1,natoms)=xt+delx(1) 189 rv(2,natoms)=yt+delx(2) 190 rv(3,natoms)=zt+delx(3) 191 if (rv(1,natoms) .lt. perlb(1)) 192 * rv(1,natoms)=rv(1,natoms)+perlen(1) 193 if (rv(1,natoms) .ge. perub(1)) 194 * rv(1,natoms)=rv(1,natoms)-perlen(1) 195 if (rv(2,natoms) .lt. perlb(2)) 196 * rv(2,natoms)=rv(2,natoms)+perlen(2) 197 if (rv(2,natoms) .ge. perub(2)) 198 * rv(2,natoms)=rv(2,natoms)-perlen(2) 199 if (rv(3,natoms) .lt. perlb(3)) 200 * rv(3,natoms)=rv(3,natoms)+perlen(3) 201 if (rv(3,natoms) .ge. perub(3)) 202 * rv(3,natoms)=rv(3,natoms)-perlen(3) 203253 continue 204252 continue 205251 continue 206250 continue 207 goto 1 208 end 209cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 210 subroutine defvalue() 211 include 'createAtoms.h' 212 integer oldtype 213 character*8 lattype 214 common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3), 215 * periodicity(3),xbound(2),ybound(2),zbound(2),strain(3), 216 * delx(3),rcell(3),ccell(nelmax) 217 common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3), 218 * avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100) 219 dimension cent(3),plane(3),planec(3),surface(3) 220 namelist /defcard/ xmin,xmax,ymin,ymax,zmin,zmax, 221 * oldtype,newtype,prob,cent,plane,dist,surface, 222 * nramp,ndirec,dismax,theta,phi,nscrew,anglemin,anglemax, 223 * xbottom,xheight,zbottom,zheight,nrepeatsx,nrepeatsz 224 data pi/3.1415927/ 225 2261 continue 227 xmin=-1.0e-20 228 xmax=1.0e20 229 ymin=-1.0e20 230 ymax=1.0e20 231 zmin=-1.0e20 232 zmax=1.0e20 233 oldtype=0 234 dist=-1.0 235 newtype=-11 236 cent(1)=0.5*(perlb(1)+perub(1)) 237 cent(2)=0.5*(perlb(2)+perub(2)) 238 cent(3)=0.5*(perlb(3)+perub(3)) 239 surface(1)=0.0 240 surface(2)=0.0 241 surface(3)=0.0 242 nramp=1 243 nscrew=0 244 anglemin=0.0 245 anglemax=360.0 246 ndirec=3 247 dismax=0.0 248 theta=5000.0 249 phi=5000.0 250 nrepeatsx=0 251 nrepeatsz=0 252 read(5,defcard) 253 if (newtype .lt. -10 .and. dismax .eq. 0.0) return 254 if (newtype .gt. ntypes) then 255 ntypes=ntypes+1 256 if (oldtype .ne. 0) then 257 amass(ntypes)=amass(oldtype) 258 else 259 amass(ntypes)=amass(ntypes-1) 260 endif 261 endif 262 if (surface(1) .eq. 1.0) then 263 perlb(1)=perlb(1)-10.0 264 perub(1)=perub(1)+10.0 265 perlen(1)=perub(1)-perlb(1) 266 endif 267 if (surface(2) .eq. 1.0) then 268 perlb(2)=perlb(2)-10.0 269 perub(2)=perub(2)+10.0 270 perlen(2)=perub(2)-perlb(2) 271 endif 272 if (surface(3) .eq. 1.0) then 273 perlb(3)=perlb(3)-10.0 274 perub(3)=perub(3)+10.0 275 perlen(3)=perub(3)-perlb(3) 276 endif 277 vol=alat(1)*alat(2)*alat(3) 278 if (dismax .ne. 0.0) then 279 astart=anglemin 280 afinish=anglemax 281 if (anglemax .lt. anglemin) then 282 anglemin=afinish 283 anglemax=astart 284 endif 285 if (nramp .eq. 1) then 286 pstart=xmin 287 pfinish=xmax 288 if (xmax .lt. xmin) then 289 xmin=pfinish 290 xmax=pstart 291 endif 292 endif 293 if (nramp .eq. 2) then 294 pstart=ymin 295 pfinish=ymax 296 if (ymax .lt. ymin) then 297 ymin=pfinish 298 ymax=pstart 299 endif 300 endif 301 if (nramp .eq. 3) then 302 pstart=zmin 303 pfinish=zmax 304 if (zmax .lt. zmin) then 305 zmin=pfinish 306 zmax=pstart 307 endif 308 endif 309 do 33 i=1,natoms 310 if (rv(1,i) .lt. xmin) goto 33 311 if (rv(1,i) .gt. xmax) goto 33 312 if (rv(2,i) .lt. ymin) goto 33 313 if (rv(2,i) .gt. ymax) goto 33 314 if (rv(3,i) .lt. zmin) goto 33 315 if (rv(3,i) .gt. zmax) goto 33 316 if (nscrew .eq. 0) then 317 if (nramp .eq. 0) then 318 rv(ndirec,i)=rv(ndirec,i)+dismax 319 else 320 rv(ndirec,i)=rv(ndirec,i)+ 321 * (rv(nramp,i)-pstart)/(pfinish-pstart)*dismax 322 endif 323 else 324 if (ndirec .eq. 1) then 325 dx=rv(2,i)-cent(2) 326 dy=rv(3,i)-cent(3) 327 dxmax=abs(ymax-cent(2)) 328 dxmin=abs(ymin-cent(2)) 329 rmax=min(dxmin,dxmax) 330 endif 331 if (ndirec .eq. 2) then 332 dx=rv(3,i)-cent(3) 333 dy=rv(1,i)-cent(1) 334 dxmax=abs(zmax-cent(3)) 335 dxmin=abs(zmin-cent(3)) 336 rmax=min(dxmin,dxmax) 337 endif 338 if (ndirec .eq. 3) then 339 dx=rv(1,i)-cent(1) 340 dy=rv(2,i)-cent(2) 341 dxmax=abs(xmax-cent(1)) 342 dxmin=abs(xmin-cent(1)) 343 rmax=min(dxmin,dxmax) 344 endif 345 dr=sqrt(dx**2+dy**2) 346 if (dr .eq. 0.0) then 347 angle=0.0 348 else 349 angle=acos(dx/dr)*180.0/pi 350 if (dy .lt. 0.0) angle=360.0-angle 351 endif 352 if (angle .lt. anglemin) goto 33 353 if (angle .gt. anglemax) goto 33 354 if (nscrew .eq. 1) then 355 astart1=astart 356 afinish1=afinish 357 else if (nscrew .eq. 2) then 358 if (astart .eq. 0.0 .and. 359 * afinish .eq. 180.0) then 360 if (dr .le. dxmax) then 361 astart1=astart 362 else 363 astart1=acos(dxmax/dr)*180.0/pi 364 endif 365 if (dr .le. dxmin) then 366 afinish1=afinish 367 else 368 afinish1=180.0-acos(dxmin/dr)*180.0/pi 369 endif 370 else if (astart .eq. 180.0 .and. 371 * afinish .eq. 0.0) then 372 if (dr .le. dxmin) then 373 astart1=astart 374 else 375 astart1=180.0-acos(dxmin/dr)*180.0/pi 376 endif 377 if (dr .le. dxmax) then 378 afinish1=afinish 379 else 380 afinish1=acos(dxmax/dr)*180.0/pi 381 endif 382 else if (astart .eq. 180.0 .and. 383 * afinish .eq. 360.0) then 384 if (dr .le. dxmin) then 385 astart1=astart 386 else 387 astart1=180.0+acos(dxmin/dr)*180.0/pi 388 endif 389 if (dr .le. dxmax) then 390 afinish1=afinish 391 else 392 afinish1=360.0-acos(dxmax/dr)*180.0/pi 393 endif 394 else if (astart .eq. 360.0 .and. 395 * afinish .eq. 180.0) then 396 if (dr .le. dxmax) then 397 astart1=astart 398 else 399 astart1=360.0-acos(dxmax/dr)*180.0/pi 400 endif 401 if (dr .le. dxmin) then 402 afinish1=afinish 403 else 404 afinish1=180.0+acos(dxmin/dr)*180.0/pi 405 endif 406 else 407 write(6,*)'problem in nscrew=2' 408 stop 409 endif 410 endif 411 if (dr .gt. rmax) dr=rmax 412 rv(ndirec,i)=rv(ndirec,i)+dr/rmax* 413 * (angle-astart1)/(afinish1-astart1)*dismax 414 endif 41533 continue 416 else if (theta .eq. 5000.0 .and. phi .eq. 5000.0 .and. 417 * nrepeatsx*nrepeatsz .eq. 0.0) then 418 if (dist .lt. 0.0) then 419 do 2 i=1,natoms 420 if (rv(1,i) .lt. xmin) goto 2 421 if (rv(1,i) .gt. xmax) goto 2 422 if (rv(2,i) .lt. ymin) goto 2 423 if (rv(2,i) .gt. ymax) goto 2 424 if (rv(3,i) .lt. zmin) goto 2 425 if (rv(3,i) .gt. zmax) goto 2 426 if (oldtype .eq. itype(i) .or. oldtype .eq. 0) then 427 if (ranl() .lt. prob) itype(i)=newtype 428 endif 4292 continue 430 else 431 if (theta .eq. 5000.0) then 432 pnorm=sqrt((plane(1)*alat(2)*alat(3))**2+ 433 * (plane(2)*alat(1)*alat(3))**2+ 434 * (plane(3)*alat(1)*alat(2))**2) 435 do 30 n=1,3 436 plane(n)=vol*plane(n)/(alat(n)*pnorm) 43730 continue 438 do 20 n=1,3 439 planec(n)=0.0 440 do 20 m=1,3 441 planec(n)=planec(n)+roter(n,m)*plane(m) 44220 continue 443 else 444 planec(1)=sin(pi/180.0*theta)*cos(pi/180.0*phi) 445 planec(3)=sin(pi/180.0*theta)*sin(pi/180.0*phi) 446 planec(2)=cos(pi/180.0*theta) 447 endif 448 do 22 i=1,natoms 449 dis1=rv(1,i)-cent(1) 450 dis2=rv(2,i)-cent(2) 451 dis3=rv(3,i)-cent(3) 452 if (dis1*planec(1)+dis2*planec(2)+dis3*planec(3) 453 * .gt. dist) then 454 if (oldtype .eq. itype(i) .or. oldtype .eq. 0) 455 * itype(i)=newtype 456 endif 45722 continue 458 endif 459 else if (nrepeatsx*nrepeatsz .ne. 0) then 460 do 4 i=1,natoms 461 if (rv(1,i) .le. xmin .or. rv(1,i) .ge. xmax) then 462 ylocalx = xbottom 463 else 464 ylocalx = xbottom + 0.5*xheight*(1.0+ 465 * cos(2.0*pi*nrepeatsx/(xmax-xmin)* 466 * (rv(1,i)-(xmin+xmax)/(2.0*nrepeatsx)))) 467 endif 468 if (rv(3,i) .le. zmin .or. rv(3,i) .ge. zmax) then 469 ylocalz = zbottom 470 else 471 ylocalz = zbottom + 0.5*zheight*(1.0+ 472 * cos(2.0*pi*nrepeatsz/(zmax-zmin)* 473 * (rv(3,i)-(zmin+zmax)/(2.0*nrepeatsz)))) 474 endif 475 ylocal = ylocalx 476 if (ylocal .gt. ylocalz) ylocal = ylocalz 477 if (rv(2,i) .gt. ylocal) then 478 itype(i)=newtype 479 endif 4804 continue 481 endif 482 483 goto 1 484 end 485cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 486 subroutine storelat(latty,avec,bvec,cvec) 487 character*8 lstored(10),latty 488 dimension avecs(3,10),bvecs(3,10),cvecs(3,10) 489 dimension avec(3),bvec(3),cvec(3) 490 data lstored/'fcc ','bcc ','sc ',7*' '/ 491 data avecs/1.0,1.0,0.0, 1.0,1.0,-1.0, 2.0,0.0,0.0, 21*0.0/ 492 data bvecs/0.0,1.0,1.0, -1.0,1.0,1.0, 0.0,2.0,0.0, 21*0.0/ 493 data cvecs/1.0,0.0,1.0, 1.0,-1.0,1.0, 0.0,0.0,2.0, 21*0.0/ 494 imatch=0 495 do 10 i=1,10 496 if (latty .ne. lstored(i)) goto 10 497 imatch = 1 498 do 5 j=1,3 499 avec(j)=avecs(j,i) 500 bvec(j)=bvecs(j,i) 501 cvec(j)=cvecs(j,i) 502 5 continue 50310 continue 504 if (imatch .eq. 1) return 505 write(6,15) latty 50615 format('could not find lattice type ') 507 stop 508 end 509cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 510 subroutine applystrain() 511 include 'createAtoms.h' 512 character*8 lattype 513 common /latinfo/ lattype,alat(3),xrot(3),yrot(3),zrot(3), 514 * periodicity(3),xbound(2),ybound(2),zbound(2),strain(3), 515 * delx(3),rcell(3),ccell(nelmax) 516 common /tmpvar/ avec(3),bvec(3),cvec(3),xold(3),yold(3),zold(3), 517 * avecp(3),bvecp(3),cvecp(3),rcellp(3),roter(3,3),rvs(3,100) 518 do i=1,3 519 periodicity(i)=(1.0+strain(i))*periodicity(i) 520 avec(i)=(1.0+strain(i))*avec(i) 521 bvec(i)=(1.0+strain(i))*bvec(i) 522 cvec(i)=(1.0+strain(i))*cvec(i) 523 enddo 524 return 525 end 526cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 527 subroutine hitvalue() 528 include 'createAtoms.h' 529 namelist /hitcard/ xmin,xmax,ymin,ymax,zmin,zmax,nhit 530 nhitcards=0 5311001 nhit=0 532 xmin=-10000.0 533 xmax=10000.0 534 ymin=-10000.0 535 ymax=10000.0 536 zmin=-10000.0 537 zmax=10000.0 538 read(5,hitcard) 539 if (nhit .eq. 0) then 540 return 541 endif 542 nchange=0 543 do 1002 i=1,natoms 544 if (rv(1,i) .lt. xmin) goto 1002 545 if (rv(1,i) .gt. xmax) goto 1002 546 if (rv(2,i) .lt. ymin) goto 1002 547 if (rv(2,i) .gt. ymax) goto 1002 548 if (rv(3,i) .lt. zmin) goto 1002 549 if (rv(3,i) .gt. zmax) goto 1002 550 nchange=nchange+1 551 nhittag(i)=nhit 5521002 continue 553 if (nchange .ne. 0) nhitcards=nhitcards+1 554 goto 1001 555 end 556cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 557 subroutine systemshift() 558 include 'createAtoms.h' 559 namelist /shiftcard/ mode 560 mode=0 561 read(5,shiftcard) 562 if (mode .eq. 0) return 563 if (mode .eq. 1) then 564 top=perub(2) 565 bottom=perlb(2) 566 sdely=0.0 567 do i=1,natoms 568 if (nhittag(i) .eq. 3 .and. top .gt. rv(2,i)) top=rv(2,i) 569 if (nhittag(i) .eq. 4 .and. bottom .lt. rv(2,i)) 570 * bottom=rv(2,i) 571 enddo 572 sdely=-0.5*(top+bottom) 573 do i=1,natoms 574 rv(2,i)=rv(2,i)+sdely 575 enddo 576 endif 577 if (mode .eq. 2) then 578c This shift only assumes a maximum of two plane spacings. 579 xmin=10.0 580 do i=1,natoms 581 if (rv(1,i) .lt. xmin) xmin=rv(1,i) 582 enddo 583 shiftneg=-10.0 584 shiftpos=10.0 585 do i=1,natoms 586 shift=rv(1,i)-xmin 587 shift=shift-perlen(1)*nint(shift/perlen(1)) 588 if (shift .gt. 0.01 .and. shift .lt. shiftpos) 589 * shiftpos=shift 590 if (shift .lt. -0.01 .and. shift .gt. shiftneg) 591 * shiftneg=shift 592 enddo 593 if (-shiftneg .gt. shiftpos) then 594 shift=0.5*shiftneg 595 else 596 shift=0.5*shiftpos 597 endif 598 perlb(1)=xmin+shift 599 perub(1)=perlb(1)+perlen(1) 600 endif 601 end 602cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 603 subroutine colect() 604 include 'createAtoms.h' 605 do i=1,natoms 606 if (rv(1,i) .lt. perlb(1)) rv(1,i)=rv(1,i)+perlen(1) 607 if (rv(1,i) .ge. perub(1)) rv(1,i)=rv(1,i)-perlen(1) 608 if (rv(2,i) .lt. perlb(2)) rv(2,i)=rv(2,i)+perlen(2) 609 if (rv(2,i) .ge. perub(2)) rv(2,i)=rv(2,i)-perlen(2) 610 if (rv(3,i) .lt. perlb(3)) rv(3,i)=rv(3,i)+perlen(3) 611 if (rv(3,i) .ge. perub(3)) rv(3,i)=rv(3,i)-perlen(3) 612 enddo 613 return 614 end 615cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 616 real function ranl() 617 common /iran/ iseed 618 iseed=iseed*125 619 iseed=iseed-2796203*(iseed/2796203) 620 ranl=abs(iseed/2796203.0) 621 return 622 end 623cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 624 subroutine delete() 625 include 'createAtoms.h' 626 nw_del=0 627 do i=1,ntypes 628 nntype(i)=0 629 enddo 630 natoms0=natoms 631 natoms=0 632 do 1 i=1,natoms0 633 if (itype(i) .eq. -1) nw_del=nw_del+1 634 if (itype(i) .le. 0) goto 1 635 natoms=natoms+1 636 ntag(natoms)=natoms 637 itype(natoms)=itype(i) 638 nhittag(natoms)=nhittag(i) 639 nntype(itype(natoms))=nntype(itype(natoms))+1 640 rv(1,natoms)=rv(1,i) 641 rv(2,natoms)=rv(2,i) 642 rv(3,natoms)=rv(3,i) 6431 continue 644 return 645 end 646cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 647 subroutine randomize() 648 include 'createAtoms.h' 649 dimension rv0(6,natmax),itype0(natmax),nhittag0(natmax) 650 do 1 i=1,natoms 651 itype0(i)=itype(i) 652 if (nhitcards .gt. 0) nhittag0(i)=nhittag(i) 653 rv0(1,i)=rv(1,i) 654 rv0(2,i)=rv(2,i) 655 rv0(3,i)=rv(3,i) 656 rv0(4,i)=rv(4,i) 657 rv0(5,i)=rv(5,i) 658 rv0(6,i)=rv(6,i) 6591 continue 660 im=1 66110 im=im*2 662 if (im .lt. natoms) goto 10 663 ia=365 664 ic=8161 665 ivalue=mod(ilatseed,natoms) 666 do 2 i=1,natoms 66720 ivalue=mod(ia*ivalue+ic,im) 668 if (ivalue .ge. natoms) goto 20 669 ntmp=ivalue+1 670 if (ntmp .ge. 1 .and. ntmp .le. natoms) then 671 ntag(ntmp)=i 672 itype(ntmp)=itype0(i) 673 if (nhitcards .gt. 0) nhittag(ntmp)=nhittag0(i) 674 rv(1,ntmp)=rv0(1,i) 675 rv(2,ntmp)=rv0(2,i) 676 rv(3,ntmp)=rv0(3,i) 677 rv(4,ntmp)=rv0(4,i) 678 rv(5,ntmp)=rv0(5,i) 679 rv(6,ntmp)=rv0(6,i) 680 endif 6812 continue 682 return 683 end 684cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 685 subroutine write0(dumpfile) 686 include 'createAtoms.h' 687 character*32 dumpfile, typ 688 open (unit=2,file=dumpfile) 689 write(2,*) natoms 690 write(2,100) 'GaN lattice' 691 do i=1,natoms 692 if (itype(i).eq.1) then 693 typ = 'Ga' 694 else 695 typ = 'N' 696 endif 697 write(2,150) typ(1:length(typ)),rv(1,i),rv(2,i),rv(3,i) 698 enddo 699 50 format(i8) 700 100 format(a) 701 150 format(a3,3f13.5) 702 close(2) 703 return 704 end 705cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 706 subroutine write1(dumpfile) 707 include 'createAtoms.h' 708 character*32 dumpfile 709 zero=0.0 710 open (unit=2,file=dumpfile) 711 write(2,*) 'ITEM: NUMBER OF ATOMS' 712 write(2,*) natoms 713 write(2,*) 'ITEM: BOX BOUNDS' 714 write(2,*) perlb(1),perub(1) 715 write(2,*) perlb(2),perub(2) 716 write(2,*) perlb(3),perub(3) 717 write(2,*) 'ITEM: TIME' 718 write(2,*) 0.0 719 write(2,*) 'ITEM: ATOMS' 720 if (nhitcards .eq. 0) then 721 do i=1,natoms 722 write(2,*) ntag(i),itype(i),rv(1,i),rv(2,i),rv(3,i) 723 enddo 724 else 725 do i=1,natoms 726 write(2,*) ntag(i),itype(i),nhittag(i), 727 * rv(1,i),rv(2,i),rv(3,i) 728 enddo 729 endif 730 close(2) 731 open(unit=2,file=dumpfile(1:index(dumpfile,' ')-1)//'.vel') 732 write(2,*) 'ITEM: NUMBER OF ATOMS' 733 write(2,*) natoms 734 write(2,*) 'ITEM: BOUNDARY VELOCITY' 735 write(2,*) zero,zero,zero 736 write(2,*) 'ITEM: TIME' 737 write(2,*) zero 738 write(2,*) 'ITEM: VELOCITIES' 739 do i=1,natoms 740 write(2,*) ntag(i),rv(4,i),rv(5,i),rv(6,i) 741 enddo 742 close(2) 743 return 744 end 745cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 746 subroutine write2(dumpfile) 747 include 'createAtoms.h' 748 data conmas/1.0365e-4/ 749 character*32 dumpfile 750 open (unit=2,file=dumpfile) 751 zero=0.0 752 write(2,*) 'using dynamo' 753 write(2,9502) natoms,ntypes,zero 7549502 format(2i10,e15.8) 755 write(2,9503) (perub(i),i=1,3),(perlb(i),i=1,3) 7569503 format(3e25.16) 757 write(2,9504) (amass(i)*conmas,ielement(i),i=1,ntypes) 7589504 format(e25.16,i10) 759 write(2,9505) ((rv(i,j),i=1,6),itype(j),j=1,natoms) 7609505 format(3e25.16/3e25.16/i10) 761 write(2,9503) zero,zero,zero 762 close(2) 763 return 764 end 765cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 766 subroutine write3(dumpfile) 767 include 'createAtoms.h' 768 character*32 dumpfile 769 open (unit=2,file=dumpfile) 770 area=(1.0-float(nw_del)/float(natoms0))*perlen(2)*perlen(3) 771 write(2,*) 'Start File for LAMMPS ',area 772 write(2,*) 773 write(2,*) natoms,' atoms' 774 write(2,*) 775 write(2,*) ntypes,' atom types' 776 write(2,*) 777 write(2,*) perlb(1),perub(1),' xlo xhi' 778 write(2,*) perlb(2),perub(2),' ylo yhi' 779 write(2,*) perlb(3),perub(3),' zlo zhi' 780 if (xy .ge. 1.0e8 .or. xz .ge. 1.0e8 .or. yz .ge. 1.0e8) then 781 write(2,*) 782 else 783 write(2,*) 784 write(2,*)xy,xz,yz,' xy xz yz' 785 write(2,*) 786 endif 787 write(2,*) 'Masses' 788 write(2,*) 789 do i=1,ntypes 790 write(2,*)i,amass(i) 791 enddo 792 write(2,*) 793 write(2,*) 'Atoms' 794 write(2,*) 795 do i=1,natoms 796 write(2,*) i,itype(i),rv(1,i),rv(2,i),rv(3,i) 797 enddo 798 write(2,*) 799 write(2,*) 'Velocities' 800 write(2,*) 801 zero=0.0 802 do i=1,natoms 803 write(2,*) i,rv(4,i),rv(5,i),rv(6,i) 804 enddo 805 close(2) 806 return 807 end 808cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 809 subroutine write4(dumpfile) 810 include 'createAtoms.h' 811 character*32 dumpfile,full 812 dumpfile='ens' 813 full=dumpfile(1:index(dumpfile,' ')-1)//'.case' 814 open(unit=2,file=full) 815 write(2,101) 816101 format('FORMAT') 817 write(2,102) 818102 format('type: ensight') 819 write(2,103) 820103 format('GEOMETRY') 821 write(2,104)'ens.case.****.geo' 822104 format('model: 1 ',a30) 823 write(2,105) 824105 format('VARIABLE') 825 write(2,111) 826111 format('TIME') 827 write(2,112) 828112 format('time set: 1') 829 write(2,113) 830113 format('number of steps: 1') 831 write(2,114) 832114 format('filename start number: 1') 833 write(2,115) 834115 format('filename increment: 1') 835 write(2,116) 836116 format('time values:') 837 write(2,117) 838117 format('1') 839 close(2) 840 full=dumpfile(1:index(dumpfile,' ')-1)//'.case.0001.geo' 841 open(unit=2,file=full) 842 write(2,201) 843 write(2,202) 844 write(2,203) 845 write(2,204) 846 write(2,205) 847 write(2,206)natoms 848 do j=1,natoms 849 write(2,207)j,rv(1,j),rv(2,j),rv(3,j) 850 enddo 851 do j=1,ntypes 852 write(2,208)j 853 write(2,209)j 854 write(2,210) 855 write(2,206)nntype(j) 856 do k=1,natoms 857 if (itype(k) .eq. j) write(2,211)k,k 858 enddo 859 enddo 860201 format('spparks input') 861202 format('geometry for element group 1') 862203 format('node id given') 863204 format('element id given') 864205 format('coordinates') 865206 format(i8) 866207 format(i8,3e12.5) 867208 format('part',i2) 868209 format('atombox',i2) 869210 format('point') 870211 format(2i8) 871 close(2) 872 return 873 end 874cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 875 subroutine write5(dumpfile) 876 include 'createAtoms.h' 877 character*32 dumpfile,full 878 character*200 sentence 879 full=dumpfile(1:index(dumpfile,' ')-1)//'.lat' 880 open (unit=2,file=full) 881 call neigen() 882 nmax=numneigh(1) 883 do i=1,natoms 884 if (nmax .lt. numneigh(i)) nmax=numneigh(i) 885 enddo 886 write(2,*) 'spparks input lattice' 887 write(2,*) 888 write(2,*) '3 dimension' 889 write(2,*) natoms,' vertices' 890 write(2,*) nmax,' max connectivity' 891 write(2,*) perlb(1),perub(1),' xlo xhi' 892 write(2,*) perlb(2),perub(2),' ylo yhi' 893 write(2,*) perlb(3),perub(3),' zlo zhi' 894 write(2,*) 895 write(2,*) 'Vertices' 896 write(2,*) 897 do i=1,natoms 898 write(2,*) i,rv(1,i),rv(2,i),rv(3,i) 899 enddo 900 write(2,*) 901 write(2,*) 'Edges' 902 write(2,*) 903 do i=1,natoms 904 write(sentence,*) i,(neigh(k,i),k=1,numneigh(i)) 905 write(2,*)sentence(1:len_trim(sentence)) 906 enddo 907 close(2) 908 full=dumpfile(1:index(dumpfile,' ')-1)//'.conf' 909 open (unit=2,file=full) 910 nvalues=1 911 write(2,*) 'spparks input configuration' 912 write(2,*) natoms,nvalues 913 write(2,*) 914 do i=1,natoms 915 write(2,*) i,itype(i) 916 enddo 917 close(2) 918 return 919 end 920cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 921 subroutine disturb() 922 include 'createAtoms.h' 923 dimension strain(3) 924 namelist /disturbcard/ dismax,strain 925 dismax=0.0 926 strain(1)=0.0 927 strain(2)=0.0 928 strain(3)=0.0 929 read(5,disturbcard) 930 if (dismax .ne. 0.0) then 931 do 1 i=1,natoms 932 rv(1,i)=rv(1,i)+dismax*(ranl()-0.5) 933 rv(2,i)=rv(2,i)+dismax*(ranl()-0.5) 934 rv(3,i)=rv(3,i)+dismax*(ranl()-0.5) 9351 continue 936 endif 937 if (strain(1)**2+strain(2)**2+strain(3)**2 .ne. 0.0) then 938 do 2 i=1,3 939 perlen(i)=perlen(i)*(1.0+strain(i)) 940 perub(i)=perlb(i)+perlen(i) 9412 continue 942 do 3 i=1,natoms 943 do 3 j=1,3 944 rv(j,i)=perlb(j)+(rv(j,i)-perlb(j))*(1.0+strain(j)) 9453 continue 946 endif 947 return 948 end 949cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 950 subroutine outputfiles() 951 character*32 dynamo,paradyn,lammps,xyz,ensight,spparks 952 namelist /filecard/ dynamo,paradyn,lammps,xyz,ensight,spparks 953 xyz="none" 954 paradyn="none" 955 dynamo="none" 956 lammps="none" 957 ensight="none" 958 spparks="none" 959 read(5,filecard) 960 if (xyz .ne. "none") call write0(xyz) 961 if (paradyn .ne. "none") call write1(paradyn) 962 if (dynamo .ne. "none") call write2(dynamo) 963 if (lammps .ne. "none") call write3(lammps) 964 if (ensight .ne. "none") call write4(ensight) 965 if (spparks .ne. "none") call write5(spparks) 966 return 967 end 968cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 969 integer function length(str) 970 character*(*) str 971 n = len(str) 972 10 if (n.gt.0) then 973 if (str(n:n).eq.' ') then 974 n = n - 1 975 goto 10 976 endif 977 endif 978 length = n 979 return 980 end 981cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 982 subroutine velgen() 983 include 'createAtoms.h' 984 data boltz/8.617e-5/,conmas/1.0365e-4/ 985 namelist /velcard/ Tbnd,Tmid,naxis,nmesh,equal,ensureTave, 986 * steeper 987 dimension Tmesh(200),pmesh(200),vcm(3,200),tmass(200),ekin(200), 988 * scale(200),ncount(200) 989 Tbnd=0.0 990 Tmid=0.0 991 naxis=1 992 nmesh=100 993 equal=1.0 994 ensureTave=0.0 995 steeper=0.0 996 read(5,velcard) 997 Tbnd=equal*Tbnd 998 Tmid=equal*Tmid 999 Tave=0.5*(Tbnd+Tmid) 1000 Tbnd=Tave+(Tbnd-Tave)*(1.0+steeper) 1001 Tmid=Tave+(Tmid-Tave)*(1.0+steeper) 1002 1003 if (Tbnd+Tmid .le. 0.0) then 1004 do i=1,natoms 1005 rv(4,i)=0.0 1006 rv(5,i)=0.0 1007 rv(6,i)=0.0 1008 enddo 1009 else 1010 wmesh=perlen(naxis)/nmesh 1011 do i=1,nmesh 1012 pmesh(i)=(i-0.5)*wmesh 1013 Tmesh(i)=Tmid+(Tbnd-Tmid)*abs(pmesh(i)-0.5*perlen(naxis))/ 1014 * (0.5*perlen(naxis)) 1015 vcm(1,i)=0.0 1016 vcm(2,i)=0.0 1017 vcm(3,i)=0.0 1018 tmass(i)=0.0 1019 ekin(i)=0.0 1020 ncount(i)=0 1021 enddo 1022 do i=1,natoms 1023 index=(rv(naxis,i)-perlb(naxis))/wmesh+1 1024 it=itype(i) 1025 vnorm=sqrt(Tmesh(index)*boltz/(amass(it)*conmas)) 1026 do j=1,3 1027 rv(j+3,i)=vnorm*rgauss() 1028 vcm(j,index)=vcm(j,index)+amass(it)*conmas*rv(j+3,i) 1029 enddo 1030 tmass(index)=tmass(index)+amass(it)*conmas 1031 ncount(index)=ncount(index)+1 1032 enddo 1033 do i=1,nmesh 1034 do j=1,3 1035 vcm(j,i)=vcm(j,i)/tmass(i) 1036 enddo 1037 enddo 1038 do i=1,natoms 1039 index=(rv(naxis,i)-perlb(naxis))/wmesh+1 1040 do j=1,3 1041 rv(j+3,i)=rv(j+3,i)-vcm(j,index) 1042 enddo 1043 enddo 1044 do i=1,natoms 1045 index=(rv(naxis,i)-perlb(naxis))/wmesh+1 1046 it=itype(i) 1047 ekin(index)=ekin(index)+0.5*amass(it)*conmas* 1048 * (rv(4,i)**2+rv(5,i)**2+rv(6,i)**2) 1049 enddo 1050 do i=1,nmesh 1051 treal=2.0/(3.0*boltz)*ekin(i)/ncount(i) 1052 if (treal .eq. 0.0) then 1053 scale(i)=0.0 1054 else 1055 scale(i)=sqrt(Tmesh(i)/treal) 1056 endif 1057 enddo 1058 do i=1,natoms 1059 index=(rv(naxis,i)-perlb(naxis))/wmesh+1 1060 do j=4,6 1061 rv(j,i)=rv(j,i)*scale(index) 1062 enddo 1063 enddo 1064 if (ensureTave .eq. 1.0) then 1065 ekint=0.0 1066 do i=1,natoms 1067 it=itype(i) 1068 ekint=ekint+0.5*amass(it)*conmas* 1069 * (rv(4,i)**2+rv(5,i)**2+rv(6,i)**2) 1070 enddo 1071 trealt=2.0/(3.0*boltz)*ekint/natoms 1072 if (trealt .eq. 0.0) then 1073 scalet=0.0 1074 else 1075 scalet=sqrt(Tave/trealt) 1076 endif 1077 do i=1,natoms 1078 do j=4,6 1079 rv(j,i)=rv(j,i)*scalet 1080 enddo 1081 enddo 1082 endif 1083 endif 1084 return 1085 end 1086cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1087 function rgauss() 1088 data twopi/6.283185308/ 10891 continue 1090 a1 = ranl() 1091 if (a1 .eq. 0.0) goto 1 1092 a2 = ranl() 1093 rgauss = sqrt(-2.0*log(a1))*cos(twopi*a2) 1094 return 1095 end 1096cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1097 subroutine neigen() 1098 include 'createAtoms.h' 1099 dimension idbox(100,100,100,100),nbox(100,100,100) 1100 cutoff=100.0 1101 xmin=rv(1,1) 1102 xmax=rv(1,1) 1103 ymin=rv(2,1) 1104 ymax=rv(2,1) 1105 zmin=rv(3,1) 1106 zmax=rv(3,1) 1107 do i=2,natoms 1108 if (xmin .gt. rv(1,i)) xmin=rv(1,i) 1109 if (xmax .lt. rv(1,i)) xmax=rv(1,i) 1110 if (ymin .gt. rv(2,i)) ymin=rv(2,i) 1111 if (ymax .lt. rv(2,i)) ymax=rv(2,i) 1112 if (zmin .gt. rv(3,i)) zmin=rv(3,i) 1113 if (zmax .lt. rv(3,i)) zmax=rv(3,i) 1114 dx=rv(1,i)-rv(1,1) 1115 dy=rv(2,i)-rv(2,1) 1116 dz=rv(3,i)-rv(3,1) 1117 dx=dx-perlen(1)*nint(dx/perlen(1)) 1118 dy=dy-perlen(2)*nint(dy/perlen(2)) 1119 dz=dz-perlen(3)*nint(dz/perlen(3)) 1120 tmp=dx**2+dy**2+dz**2 1121 if (tmp .lt. cutoff) cutoff=tmp 1122 enddo 1123 cutoff=sqrt(cutoff)+0.1 1124 nx=(xmax-xmin)/cutoff+1 1125 ny=(ymax-ymin)/cutoff+1 1126 nz=(zmax-zmin)/cutoff+1 1127 if (nx .gt. 100) nx=100 1128 if (ny .gt. 100) ny=100 1129 if (nz .gt. 100) nz=100 1130 delx=(xmax-xmin)/nx 1131 dely=(ymax-ymin)/ny 1132 delz=(zmax-zmin)/nz 1133 do n1=1,nx 1134 do n2=1,ny 1135 do n3=1,nz 1136 nbox(n1,n2,n3)=0 1137 enddo 1138 enddo 1139 enddo 1140 do j=1,natoms 1141 numneigh(j)=0 1142 n1=(rv(1,j)-xmin)/delx+1 1143 if (n1 .lt. 1) n1=1 1144 if (n1 .gt. nx) n1=nx 1145 n2=(rv(2,j)-ymin)/dely+1 1146 if (n2 .lt. 1) n2=1 1147 if (n2 .gt. ny) n2=ny 1148 n3=(rv(3,j)-zmin)/delz+1 1149 if (n3 .lt. 1) n3=1 1150 if (n3 .gt. nz) n3=nz 1151 nbox(n1,n2,n3)=nbox(n1,n2,n3)+1 1152 idbox(n1,n2,n3,nbox(n1,n2,n3))=j 1153 enddo 1154 do j=1,natoms 1155 n1=(rv(1,j)-xmin)/delx+1 1156 if (n1 .lt. 1) n1=1 1157 if (n1 .gt. nx) n1=nx 1158 n2=(rv(2,j)-ymin)/dely+1 1159 if (n2 .lt. 1) n2=1 1160 if (n2 .gt. ny) n2=ny 1161 n3=(rv(3,j)-zmin)/delz+1 1162 if (n3 .lt. 1) n3=1 1163 if (n3 .gt. nz) n3=nz 1164 do na=n1-1,n1+1 1165 if (na .lt. 1) then 1166 naa=nx 1167 else if (na .gt. nx) then 1168 naa=1 1169 else 1170 naa=na 1171 endif 1172 do nb=n2-1,n2+1 1173 if (nb .lt. 1) then 1174 nbb=ny 1175 else if (nb .gt. ny) then 1176 nbb=1 1177 else 1178 nbb=nb 1179 endif 1180 do nc=n3-1,n3+1 1181 if (nc .lt. 1) then 1182 ncc=nz 1183 else if (nc .gt. nz) then 1184 ncc=1 1185 else 1186 ncc=nc 1187 endif 1188 do 100 ndd=1,nbox(naa,nbb,ncc) 1189 k=idbox(naa,nbb,ncc,ndd) 1190 if (k .ge. j) goto 100 1191 dx=rv(1,k)-rv(1,j) 1192 dy=rv(2,k)-rv(2,j) 1193 dz=rv(3,k)-rv(3,j) 1194 dx=dx-perlen(1)*nint(dx/perlen(1)) 1195 dy=dy-perlen(2)*nint(dy/perlen(2)) 1196 dz=dz-perlen(3)*nint(dz/perlen(3)) 1197 tmp=sqrt(dx**2+dy**2+dz**2) 1198 if (tmp .le. cutoff) then 1199 numneigh(j)=numneigh(j)+1 1200 neigh(numneigh(j),j)=k 1201 numneigh(k)=numneigh(k)+1 1202 neigh(numneigh(k),k)=j 1203 endif 1204100 continue 1205 enddo 1206 enddo 1207 enddo 1208 enddo 1209 return 1210 end 1211