program twist
c     Create a symmetric twist GB from a perfect crystal.
c     Detlef Conrad 07/14/97.
c
c     General: read coordinates of a structure (make sure that
c     first N/2 atoms are in the lower part), 
c     check whether a desired Sigma is possible,
c     enlarge the lower part of the structure and rotate it 
c     with respect to the old coordinate system,
c     enlarge the upper part and rotate it with respect to the
c     lower part.
c     For some values of sigma, the code is not very stable. There
c     are atoms missing or too much atoms which is probably due to
c     truncation errors. However, up to now I found it more easy
c     to look at the improper structures and fix it by hand than
c     to make the program better.
c     units are in Angstrom.
 
      implicit real*8(a-h,o-z)
      character text*36,ltype*2,lt1*2,pdbfile*11   
      common/para/xl,yl,x0,y0,z1,z2
      common/xyz/ pos(3,80),ltype(80),zl
      dimension pos1(3,30000),pos2(3,300000)    
      dimension lt1(30000)
      character atom*6,atm*2,at*1

      pdbfile='test.pdb'
      call pdbin(pdbfile,nmol,pos,ltype,zl)     ! read structure
      xl=1.0
      yl=1.0
      a0=4.32                                   ! lattice constant of SiC
                                                ! must be changed for other structures   
      ianz=nmol                                 ! number of atoms
   
      pi=4.0*datan(1.0d0)
      text=' 90.000 90.000 90.000 P 1           1'
      open(2,file='pos1.pdb')
      open(3,file='pos2.pdb')
112   format('ATOM  ',i5,1x,a2,11x,a1,4x,3f8.3)     
113   format(a6,3(1X,f8.3),a37) 
114   format(a6,i5,1x,a2,11x,a1,4x,3f8.3)

    
      write(*,*) 'N='                            ! read desired Sigma
      read(*,*) n
      xn=n                                       ! use real*8 
      xn2=xn*xn
c     Check whether the condition Sigma*Sigma = N1*N1 + N2*N2
c     can be fulfilled. This condition is necessary and sufficient
c     for a GB with the desired Sigma for a square unit cell
c     (cf. the book of Bollmann).
c     I need real*8, that's why x0=N1, y0=N2.
      x0=0.0
      do 100,i=1,n-1                             
         x0=x0+1.0                                  
         y0=0.0
         do 110,j=1,i
            y0=y0+1.0
            z2=x0*x0+y0*y0
            if (dabs(z2-xn2).lt.0.001) goto 120
110      continue
100   continue
c     If the condition can be fulfilled, find out the values of N1 (x0)
c     and N2 (y0). I do not remember why I'm doing this again since I
c     could use the values of the previous loop, but I do not like to
c     change it.
120   if (dabs(z2-xn2).lt.0.001) then
          theta=dacos(x0/xn)
          write(*,*) 'Theta=',theta*180.0/pi   ! This is the angle of misorientation
          x0=0.0
          do 130,i=1,n-1
            x0=x0+1.0
            y0=0.0
            do 140,j=1,i
              y0=y0+1.0
              z=x0*x0+y0*y0
              if (dabs(z-xn).lt.0.0001) goto 150
140      continue
130   continue
150   if (dabs(z-xn).lt.0.001) write(*,*) 'x0,y0,z',x0,y0,z !Check: x0*x0+y0*y0=z=xn2
        g1=dasin(y0/dsqrt(xn))
        g2=dacos(x0/dsqrt(xn))
        write(*,*) 'gamma=',g2*180.0/pi
        z1=dsqrt(xn)/dsin(g2)
        z2=dsqrt(xn)/dcos(g2)
        write(*,*) 'z1,z2',z1,z2
c       Make the lower part of the structure Sigma times larger.
c       Up to now, no rotation is needed but only atoms are
c       considered which are in a square that does not coincide
c       with the original cartesian coordinate system.
        nmol=0
        j1=-nint(y0)-1
        j2=nint(x0)+1
        i2=nint(x0+y0)+1
        write(*,*) 'j1,j2,i2',j1,j2,i2
        do 200,i=0,i2
           do 210,j=j1,j2
              do 220,k=1,ianz/2
              x=i*xl+pos(1,k)               ! coordinates of the new atom
              y=j*yl+pos(2,k)
              if (insquare(x,y).eq.1) then  ! check if this is in the square
                nmol=nmol+1
                pos1(1,nmol)=x*a0
                pos1(2,nmol)=y*a0
                pos1(3,nmol)=pos(3,k) 
                lt1(nmol)=ltype(k)          ! Type of atom.
               endif
220            continue
210        continue
200     continue
      write(*,*) 'nmol=',nmol ! Check number of atoms. Must be ianz/2 * Sigma.
      
      sinus=-y0/dsqrt(xn)
      cosin=x0/dsqrt(xn)
c     gamma is the angle of which the lower part of the structure is
c     rotated with respect to the old coordinate system,
c     this has nothing to do with the misorientation !
      gamma=dacos(cosin)
      xll=dsqrt(xn)*a0           ! new cell dimensions
      yll=xll                    ! valid only for square lattice
      zll=zl                     ! z dimension not changed 
      write(2,113) 'CRYST1 ',xll,yll,zll,text
      write(3,113) 'CRYST1 ',xll,yll,zll,text

c     This loop rotates the lower part of the new structure so that
c     the sides of the square fit into the old cartesian coordinate system.
c     The angle has nothing to do with the angle of misorientation !!! 
      do 300,i=1,nmol
         x=pos1(1,i)
         y=pos1(2,i)
         pos1(1,i)=x*cosin+y*sinus
         pos1(2,i)=-x*sinus+y*cosin
300   continue  
c     Write the lower part to file pos1.pdb. In principle, one can write
c     the whole structue in one file, but in this way it is easier to
c     check where atoms are lost (if any).
      do 310,i=1,nmol
         write(2,112) i,lt1(i),'1',pos1(1,i),pos1(2,i),pos1(3,i)
310   continue
 
c     Now enlarge the upper part of the old structure.
      nhilf=0 
      do 350,kx=-i2,i2
         do 360,ky=-i2,i2
            do 370,i=ianz/2+1,ianz
            nhilf=nhilf+1
            pos2(1,nhilf)=kx*a0+pos(1,i)*a0
            pos2(2,nhilf)=ky*a0+pos(2,i)*a0
            pos2(3,nhilf)=pos(3,i)
            lt1(nhilf) = ltype(i)
370       continue
360     continue
350   continue
c     Now rotate the upper part of the structure by an angle of
c     theta+gamma with respect to the old coordinate system. 
c     If an atom is situated in the new box, write it to file
c     pos2.pdb
      sinus=-dsin(theta+gamma)
      cosin=dcos(theta+gamma)
      nn=nmol
      do 400,i=1,nhilf
         x=pos2(1,i)
         y=pos2(2,i)
         pos2(1,i)=x*cosin+y*sinus+xl
         pos2(2,i)=-x*sinus+y*cosin+yl
         if ((pos2(1,i).le.xll).and.(pos2(1,i).ge.0.0)) then
           if ((pos2(2,i).le.yll).and.(pos2(2,i).ge.0.0)) then
           nn=nn+1
           write(3,112) nn,lt1(i),'1',pos2(1,i),pos2(2,i),pos2(3,i)
           endif
      endif
400   continue
      write(*,*) nn            ! Check: nn=Sigma * ianz/2 
      close(1)
      close(2)
      endif
c     If everything is okay, one can remove the first line from file pos2.pdb
c     and 'cat' the files pos1.pdb pos2.pdb.
      end   

      integer function insquare(x,y)
c     This is to check whether an atom is in a square which
c     sides do not coincide with the coordinate system. 
c     Coordinates of the atom is (x,y), z is not needed.
      implicit real*8(a-h,o-z)
      common/para/xl,yl,x0,y0,z1,z2
      insquare=0
      if ((y.le.(x0/y0)*x).and.(y.gt.(x0/y0)*x-z1)) then
         if ((y.ge.(-(y0/x0)*x)).and.(y.lt.(-(y0/x0)*x+z2))) insquare=1
      endif
      end   

      subroutine pdbin(pdbfile,nmol,pos,ltype,zl) 
      implicit real*8(a-h,o-z)
      character*2 ltype
      dimension pos(3,80),ltype(80)
c
      character*60 txt,pdbfile*11
      character  atom8*8,atom4*4
      character*2  atm,atmadd

      OPEN(UNIT=99,file=pdbfile,status='OLD')
      read(99,22)atom8,xl,yl,zl,txt
  22  format(a8,3f8.3,a60)
      write(*,22)atom8,xl,yl,zl,txt
      katom=20000
      nmol=0
      do26   i=1,katom
      read(99,24,err=26,end=26)atom4,j,atm,atmadd,a4,a5,a6
  24  format(a4,2x,i5,1x,a2,10x,a2,4x,3f8.3)
      if(atom4 .ne.'ATOM')goto26
       pos(1,i)=a4/xl               ! relative coordinates for x and y ...
       pos(2,i)=a5/xl
       pos(3,i)=a6                  ! ... but not for z.
       ltype(i)=atm
       
       nmol=nmol+1
  26   continue
       write(*,*) 'NMOL=',nmol
       return
       end