C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C-->   subroutine "writeeng" is used to                      %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        subroutine writeeng( dataxx,
     &                       engxy,
     &                       dataxy,
     &                       fitxy,
     &                       vora,
     &                       engerror,
     &                       dimode,
     &                       spgnum,
     &                       spgmode,
     &                       nelastic)

        use constant
        implicit none

        integer i,j
        integer dimode, spgmode, nelastic, spgnum
        integer engerror
        dimension engerror(21)

        character*40,dimension(:),allocatable:: DefMat

        double precision vora
        double precision dataxx, engxy, dataxy, fitxy
        dimension dataxx(ndef), engxy(21,ndef)
        dimension dataxy(21,ndef), fitxy(21,ndef)

        open(unit = 29, file = 'EPVDAT', status = 'unknown')

        write(29,'(2I4, f12.4)') spgnum, nelastic, vora

        allocate(DefMat(nelastic))

C       For three dimonsion structure        

        if ( dimode .EQ. 3 ) THEN

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

                if (spgmode .EQ. 1) THEN
                                           
             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'
     &                                                /),(/nelastic/))

                endif

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

                if (spgmode .EQ. 2) THEN

             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'   
     &                                                /),(/nelastic/))


                endif

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

                if (spgmode .EQ. 3) THEN

             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'
     &                                                /),(/nelastic/))


                endif

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

                if (spgmode .EQ. 4) THEN

            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'
     &                                                /),(/nelastic/))

                endif

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

                if (spgmode .EQ. 5) THEN

            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'
     &                                               /),(/nelastic/))


                endif

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

                if (spgmode .EQ. 6) THEN

            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'
     &                                               /),(/nelastic/))


                endif

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

                if (spgmode .EQ. 7) THEN

            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'
     &                                               /),(/nelastic/))


                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 (spgmode .EQ. 8) THEN

             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'
     &                                                /),(/nelastic/))


                endif

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

                if (spgmode .EQ. 9) THEN

            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'
     &                                               /),(/nelastic/))

                endif
        endif

C       For two dimonsion structure

        if ( dimode .EQ. 2 ) THEN

C       Oblique system
C       2D space group: p1, p2
C       (6 independent elastic constants)  

                if (spgmode .EQ. 1) THEN

             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, 0.0, 0.0, 0.0, strain',
     &                           'strain, strain, 0.0, 0.0, 0.0, 0.0',
     &                           'strain, 0.0, 0.0, 0.0, 0.0, strain',
     &                           '0.0, strain, 0.0, 0.0, 0.0, strain'
     &                                                /),(/nelastic/))

                endif

C       Rectanglar system
C       2D space group: pm, pg, p2mm, p2gg, p2mg, cm, c2mm      
C       (4 independent elastic constants) 

                if (spgmode .EQ. 2) THEN

             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, 0.0, 0.0, 0.0, strain',
     &                           'strain, strain, 0.0, 0.0, 0.0, 0.0'
     &                                                /),(/nelastic/))


                endif

C       Square system
C       2D space group: p4, p4mm, p4gm
C       (3 independent elastic constants)

                if (spgmode .EQ. 3) THEN

             DefMat(:)=reshape((/'strain, strain, 0.0, 0.0, 0.0, 0.0',
     &                           '0.0, 0.0, 0.0, 0.0, 0.0, strain',
     &                           'strain, 0.0, 0.0, 0.0, 0.0, 0.0'
     &                                                /),(/nelastic/))

                endif

C       Hexagon system
C       2D space group: p3, p6, p3m1, p31m, p6mm
C       (2 independent elastic constants)

                if (spgmode .EQ. 4) THEN

             DefMat(:)=reshape((/'strain, strain, 0.0, 0.0, 0.0, 0.0',
     &                           '0.0, 0.0, 0.0, 0.0, 0.0, strain'
     &                                                /),(/nelastic/))


                endif                
        endif


        do i=1,nelastic
           
           write(29,*)
           write(29,*) DefMat(i)
           write(29,*)

           if (engerror(i) .EQ. 1) THEN
                write(29,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
                write(29,*) '%%              >>> ERROR <<<           %%'
                write(29,*) '%% The structure is not at equilibrium. %%'
                write(29,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
                write(29,*)
           endif

           if (engerror(i) .EQ. 2) THEN
                write(29,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
                write(29,*) '%%             >>> ERROR <<<            %%'
                write(29,*) '%% The unstable distortion mode exists. %%'
                write(29,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
                write(29,*)
           endif

           if (engerror(i) .EQ. -1) THEN
                write(29,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
                write(29,*) '%%             >>> ERROR <<<            %%'
                write(29,*) '%%  Explained Sum of Squares > 0.1, it  %%'
                write(29,*) '%%  does not show quadratic relation.   %%'
                write(29,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
                write(29,*)
           endif

           do j=1,ndef
              write(29,'(4f15.8)') dataxx(j), engxy(i,j), 
     &                       dataxy(i,j), fitxy(i,j)
           enddo
        enddo

        deallocate(DefMat)

        close(29)

        open(unit=19,position='Append',FILE='RELAS') 
        
        do i=1,nelastic

           if (engerror(i) .EQ. 1) THEN
                write(19,*) 
                write(19,'(a4,i3,a13)') 'FOR', i, 'th distortion'
                write(19,*) '            >>> ERROR <<<             '
                write(19,*) ' The structure is not at equilibrium. '
                write(19,*)
           endif

           if (engerror(i) .EQ. 2) THEN
                write(19,*) 
                write(19,'(a4,i3,a13)') 'FOR', i, 'th distortion'
                write(19,*) '             >>> ERROR <<<            '
                write(19,*) ' The unstable distortion mode exists. '
                write(19,*)
           endif

           if (engerror(i) .EQ. -1) THEN
                write(19,*) 
                write(19,'(a4,i3,a13)') 'FOR', i, 'th distortion'
                write(19,*) '             >>> ERROR <<<            '
                write(19,*) '  Explained Sum of Squares > 0.1, it  '
                write(19,*) '  does not show quadratic relation.   '
                write(19,*)
           endif

        enddo

        close(19)

        end
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C%%%%                      End of writeeng.f                   %%%%
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
