C--------------------------------------------------------------%%%%
C-->     Subroutine "generate" is used to generate the files   %%%%
C-->     with various kinds of deformation used in the first   %%%%
C-->     principle calculation                                 %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine generate

        use constant
        implicit none             
        
        integer spgmode, nelastic
        integer spgnum, ntype

        character*10 inpos

        inpos='INPOS'
        call symmetry(spgnum, inpos)

C       Triclinic system
C       Point group classes: 1, -1
C       (21 independent elastic constants)      

        if (spgnum .GE. 1 .AND. spgnum .LE. 2) then
                nelastic=21
                spgmode=1
        endif

C       Monoclinic system
C       Point group classes: m, 2, 2/m
C       (13 independent elastic constants)  

        if (spgnum .GE. 3 .AND. spgnum .LE. 15) then
                nelastic=13
                spgmode=2
        endif

C       Orthorhombic system
C       Point group classes: 222, mm2, 2/m2/m2/m
C       (9 independent elastic constants)  

        if (spgnum .GE. 16 .AND. spgnum .LE. 74) then
                nelastic=9
                spgmode=3
        endif

C       Tetragonal II system
C       Point group classes: 4, -4, 4/m
C       (7 independent elastic constants)

        if (spgnum .GE. 75 .AND. spgnum .LE. 88) then
               nelastic=7
               spgmode=4
        endif

C       Tetragonal I system
C       Point group classes: 422, 4mm, -42m, 4/m2/m/2/m
C       (6 independent elastic constants)

        if (spgnum .GE. 89 .AND. spgnum .LE. 142) then
                nelastic=6
                spgmode=5
        endif

C       Trigonal II system
C       Point group classes: 3, -3
C       (7 independent elastic constants)         

        if (spgnum .GE. 143 .AND. spgnum .LE. 148) then
                nelastic=7
                spgmode=6
        endif

C       Trigonal I system
C       Point group classes: 32, 3m, -32/m
C       (6 independent elastic constants)   

        if (spgnum .GE. 149 .AND. spgnum .LE. 167) then
                nelastic=6
                spgmode=7
        endif

C       Hexagonal system
C       Point group classes: 6, -6, 6/m, 622, 6mm, -62m, 6/m2/m2/m
C       (5 independent elastic constants)  

        if (spgnum .GE. 168 .AND. spgnum .LE. 194) then
                nelastic=5
                spgmode=8
        endif

C       Cubic system
C       Point group classes: 23, 2/m-3, 432, -43m, 4/m-32/m
C       (3 independent elastic constants)

        if (spgnum .GE. 195 .AND. spgnum .LE. 230) then
                nelastic=3
                spgmode=9
        endif

        call recell

        call writedef(3, 
     &                spgmode, 
     &                nelastic, 
     &                spgnum)        

        select case (spgmode)
               case(1)
                 call Tricg(nelastic)
               case(2)
                 call Monocg(nelastic)
               case(3)
                 call Orthorg(nelastic)
               case(4)
                 call Tetrag7g(nelastic)
               case(5)
                 call Tetrag6g(nelastic)
               case(6)
                 call Trig7g(nelastic)
               case(7)
                 call Trig6g(nelastic)
               case(8)
                 call Hexagg(nelastic)
               case(9)
                 call Cubicg(nelastic)                
        end select

        end
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                   End of mian program                   %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Triclinic system
C%%%%%               (21 independent elastic constants)
C%%%%%                 c11  c12  c13  c14  c15  c16
C%%%%%                 c12  c22  c23  c24  c25  c26
C%%%%%                 c13  c23  c33  c34  c35  c36
C%%%%%                 c14  c24  c34  c44  c45  c46
C%%%%%                 c15  c25  c35  c45  c55  c56
C%%%%%                 c16  c26  c36  c46  c56  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
        subroutine Tricg(nelastic)

        use constant
        implicit none                
        
        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, 0.0, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, 0.0, 0.0, 
     &                           0.0, 0.0, 0.0, 0.0, strain, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           strain, 0.0, strain, 0.0, 0.0, 0.0,
     &                           strain, 0.0, 0.0, strain, 0.0, 0.0,
     &                           strain, 0.0, 0.0, 0.0, strain, 0.0,
     &                           strain, 0.0, 0.0, 0.0, 0.0, strain,
     &                           0.0, strain, strain, 0.0, 0.0, 0.0,
     &                           0.0, strain, 0.0, strain, 0.0, 0.0,
     &                           0.0, strain, 0.0, 0.0, strain, 0.0, 
     &                           0.0, strain, 0.0, 0.0, 0.0, strain,
     &                           0.0, 0.0, strain, strain, 0.0, 0.0,
     &                           0.0, 0.0, strain, 0.0, strain, 0.0,
     &                           0.0, 0.0, strain, 0.0, 0.0, strain,
     &                           0.0, 0.0, 0.0, strain, strain, 0.0,
     &                           0.0, 0.0, 0.0, strain, 0.0, strain,
     &                           0.0, 0.0, 0.0, 0.0, strain, strain
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Monoclinic system
C%%%%%               (13 independent elastic constants)
C%%%%%                 c11  c12  c13    0  c15    0
C%%%%%                 c12  c22  c23    0  c25    0
C%%%%%                 c13  c23  c33    0  c35    0
C%%%%%                   0    0    0  c44    0  c46
C%%%%%                 c15  c25  c35    0  c55    0
C%%%%%                   0    0    0  c46    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        subroutine Monocg(nelastic)

        use constant
        implicit none        

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, 0.0, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, 0.0, 0.0, 
     &                           0.0, 0.0, 0.0, 0.0, strain, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, strain, strain, 0.0, 0.0, 0.0,
     &                           strain, 0.0, strain, 0.0, 0.0, 0.0,
     &                           strain, 0.0, 0.0, 0.0, strain, 0.0,
     &                           0.0, strain, 0.0, 0.0, strain, 0.0,
     &                           0.0, 0.0, strain, 0.0, strain, 0.0,
     &                           0.0, 0.0, 0.0, strain, 0.0, strain   
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Orthorhombic system
C%%%%%               (9 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0    0
C%%%%%                 c12  c22  c23    0    0    0
C%%%%%                 c13  c23  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c55    0
C%%%%%                   0    0    0    0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        subroutine Orthorg(nelastic)

        use constant
        implicit none 

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, 0.0, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, 0.0, 0.0, 
     &                           0.0, 0.0, 0.0, 0.0, strain, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, strain, strain, 0.0, 0.0, 0.0,
     &                           strain, 0.0, strain, 0.0, 0.0, 0.0
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Tetragonal II system
C%%%%%               (7 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0  c16
C%%%%%                 c12  c11  c13    0    0 -c16
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                 c16 -c16    0    0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
        subroutine Tetrag7g(nelastic)

        use constant
        implicit none                 
        
        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, strain, 0.0,
     &                           strain, strain, strain, 0.0, 0.0, 0.0,
     &                           0.0, strain, strain, 0.0, 0.0, 0.0,
     &                           strain, 0.0, 0.0, 0.0, 0.0, strain
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Tetragonal I system
C%%%%%               (6 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0    0
C%%%%%                 c12  c11  c13    0    0    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                   0    0    0    0    0  c66
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
        subroutine Tetrag6g(nelastic)

        use constant
        implicit none                 

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, strain, 0.0,
     &                           strain, strain, strain, 0.0, 0.0, 0.0,
     &                           0.0, strain, strain, 0.0, 0.0, 0.0
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Rhombohedral II system
C%%%%%               (7 independent elastic constants)
C%%%%%                 c11  c12  c13  c14  c15    0
C%%%%%                 c12  c11  c13 -c14 -c15    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                 c14 -c14    0  c44    0 -c15
C%%%%%                 c15 -c15    0    0  c44  c14
C%%%%%                   0    0    0 -c15  c14  c66=(c11-c12)/2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine Trig7g(nelastic)

        use constant
        implicit none         

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, strain, 0.0,
     &                           strain, strain, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, strain, strain,
     &                           0.0, strain, 0.0, 0.0, 0.0, strain
     &                                                /),(/6,nelastic/))


           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                      Rhombohedral I system
C%%%%%               (6 independent elastic constants)
C%%%%%                 c11  c12  c13  c14    0    0
C%%%%%                 c12  c11  c13 -c14    0    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                 c14 -c14    0  c44    0    0
C%%%%%                   0    0    0    0  c44  c14
C%%%%%                   0    0    0    0  c14  c66=(c11-c12)/2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        subroutine Trig6g(nelastic)

        use constant
        implicit none 

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, strain, 0.0,
     &                           strain, strain, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, strain, strain
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                       Hexagonal system
C%%%%%               (5 independent elastic constants)
C%%%%%                 c11  c12  c13    0    0    0
C%%%%%                 c12  c11  c13    0    0    0
C%%%%%                 c13  c13  c33    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                   0    0    0    0    0  c66=(c11-c12)/2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine Hexagg(nelastic)

        use constant
        implicit none         

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, 0.0, 0.0, strain,
     &                           0.0, 0.0, strain, 0.0, 0.0, 0.0,
     &                           0.0, 0.0, 0.0, strain, strain, 0.0,
     &                           strain, strain, strain, 0.0, 0.0, 0.0
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                          Cubic system
C%%%%%               (3 independent elastic constants)
C%%%%%                 c11  c12  c12    0    0    0
C%%%%%                 c12  c11  c12    0    0    0
C%%%%%                 c12  c12  c11    0    0    0
C%%%%%                   0    0    0  c44    0    0
C%%%%%                   0    0    0    0  c44    0
C%%%%%                   0    0    0    0    0  c44
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        subroutine Cubicg(nelastic)

        use constant
        implicit none         

        integer i, j, k
        integer nelastic
        real DefVect, strain
        dimension DefVect(6)
        real,dimension(:,:),allocatable:: DefMat        

        do i=1,ndef
           strain=defdata(i)
           allocate(DefMat(6,nelastic))
           DefMat(:,:)=reshape((/0.0, 0.0, 0.0, strain, strain, strain,
     &                           strain, strain, 0.0, 0.0, 0.0, 0.0,
     &                           strain, strain, strain, 0.0, 0.0, 0.0
     &                                                /),(/6,nelastic/))

           do j=1,nelastic
              do k=1,6
                 DefVect(k)=DefMat(k,j)
              enddo
              call def4elas(DefVect, j, i)
           enddo
           deallocate(DefMat)
        enddo
        end
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%%                     End of generate.f                   %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

