geo.f90 8.77 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9
MODULE m_geo
  USE m_juDFT
CONTAINS
10 11
  SUBROUTINE geo(atoms,sym,cell,oneD,vacuum,input_in,tote,&
                 forcetot)
12 13 14

    !    *********************************************************************
    !    * calculates the NEW atomic positions after the results%force calculation   *
15
    !    * SUBROUTINE is based on a BFGS method implemented by M. Weinert    *
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
    !    *                                [cf. PRB 52 (9) p. 6313 (1995)]    * 
    !    *                                                                   *
    !    * as a first step we READ in the file 'inp' WITH some additional    *
    !    * information (SUBROUTINE rw_inp)                                   *
    !    * THEN recover the old geometry optimisation information from file  *
    !    * 'forces.dat' (SUBROUTINE bfsg0)                                   *
    !    * this input together WITH the NEW forces (forcetot) are now used   *
    !    * to calculate the NEW atomic positions (SUBROUTINE bfsg)           *
    !    * finally the NEW 'inp' file is written (SUBROUTINE rw_inp)         *
    !    *                                                           Gustav  *
    !
    ! input: 
    !        ntype .... total number of atom types
    !        thetad ... approx. debye temperature
    !        zat(ntype) mass number of the atom (or atomic number)
    !        xa ....... mixing factor 
    !        epsdisp .. limit for displacement to be converged
    !        epsforce . the same for force
    !        istepnow . steps to be done in this run
    !
    !    *********************************************************************
37
    USE m_types
38 39 40
    USE m_rwinp
    USE m_bfgs
    USE m_bfgs0
41
    USE m_constants
42 43
    USE m_rinpXML
    USE m_winpXML
44
    USE m_init_wannier_defaults
45 46 47 48 49
    IMPLICIT NONE
    TYPE(t_oneD),INTENT(IN)   :: oneD
    TYPE(t_cell),INTENT(IN)   :: cell
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_atoms),INTENT(IN)  :: atoms
50 51
    TYPE(t_vacuum),INTENT(IN) :: vacuum
    TYPE(t_input),INTENT(IN)  :: input_in
52 53 54 55 56
    ! ..
    ! ..  Scalar Arguments ..
    REAL,    INTENT (IN) :: tote
    ! ..
    ! ..  Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
57
    REAL,    INTENT (INOUT) :: forcetot(3,atoms%ntype)
58 59 60
    ! ..
    ! ..  Local Scalars ..
    INTEGER i,j,na ,istep0,istep,itype,jop,ieq
61
    LOGICAL lconv
62 63 64
    TYPE(t_atoms)  :: atoms_new
    ! ..
    ! ..  Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
65 66
    REAL xold(3*atoms%ntype),y(3*atoms%ntype),h(3*atoms%ntype,3*atoms%ntype),zat(atoms%ntype)
    REAL tau0(3,atoms%ntype),tau0_i(3,atoms%ntype) 
67 68 69

    TYPE(t_input):: input

70 71 72 73 74 75 76 77 78 79 80 81 82
    ! temporary variables for XML IO
    TYPE(t_input)                 :: input_temp
    TYPE(t_dimension)             :: dimension_temp
    TYPE(t_atoms)                 :: atoms_temp
    TYPE(t_cell)                  :: cell_temp
    TYPE(t_stars)                 :: stars_temp
    TYPE(t_sym)                   :: sym_temp
    TYPE(t_noco)                  :: noco_temp
    TYPE(t_vacuum)                :: vacuum_temp
    TYPE(t_sliceplot)             :: sliceplot_temp
    TYPE(t_banddos)               :: banddos_temp
    TYPE(t_obsolete)              :: obsolete_temp
    TYPE(t_enpara)                :: enpara_temp
83
    CLASS(t_xcpot),ALLOCATABLE    :: xcpot_temp
84 85 86 87
    TYPE(t_results)               :: results_temp
    TYPE(t_kpts)                  :: kpts_temp
    TYPE(t_hybrid)                :: hybrid_temp
    TYPE(t_oneD)                  :: oneD_temp
88
    TYPE(t_coreSpecInput)         :: coreSpecInput_temp
89
    TYPE(t_wann)                  :: wann_temp
90
    LOGICAL                       :: l_kpts_temp
91 92 93 94 95 96 97 98 99 100 101 102 103
    INTEGER                       :: numSpecies
    INTEGER                       :: div(3)
    INTEGER, ALLOCATABLE          :: xmlElectronStates(:,:)
    INTEGER, ALLOCATABLE          :: atomTypeSpecies(:)
    INTEGER, ALLOCATABLE          :: speciesRepAtomType(:)
    REAL, ALLOCATABLE             :: xmlCoreOccs(:,:,:)
    LOGICAL, ALLOCATABLE          :: xmlPrintCoreStates(:,:)
    CHARACTER(len=3), ALLOCATABLE :: noel_temp(:)
    CHARACTER(len=4)              :: namex_temp
    CHARACTER(len=12)             :: relcor_temp
    CHARACTER(LEN=20)             :: filename
    REAL                          :: a1_temp(3),a2_temp(3),a3_temp(3)
    REAL                          :: scale_temp, dtild_temp
104
    CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
105
    input=input_in
106 107
    atoms_new=atoms

108 109 110 111 112
    istep0 = 0
    xold = 0.0
    y = 0.0
    h= 0.0

113 114
    na = 1
    DO i = 1,atoms_new%ntype
115
       zat(i)=real(atoms%nz(i))
116
       tau0(:,i)=atoms%pos(:,na)
117 118
       na = na + atoms_new%neq(i)
    END DO
119

120
    CALL bfgs0(atoms%ntype, istep0,xold,y,h)
121 122 123

    DO itype=1,atoms%ntype
       IF (atoms%l_geo(itype)) THEN
124
          WRITE (6,'(6f10.5)') (tau0(j,itype),j=1,3), (forcetot(i,itype),i=1,3)
125 126 127 128
          DO i = 1,3
             forcetot(i,itype)=forcetot(i,itype)*REAL(atoms%relax(i,itype))
          ENDDO
          WRITE (6,'(6f10.5,a,3i2)') (tau0(j,itype),j=1,3),&
129
               &     (forcetot(i,itype),i=1,3),' atoms%relax: ', (atoms%relax(i,itype),i=1,3)
130 131 132 133 134 135 136 137
       ELSE
          DO i = 1,3
             forcetot(i,itype)=0.0
          ENDDO
       ENDIF
    ENDDO

    istep = 1
138
    CALL bfgs(atoms%ntype,istep,istep0,forcetot,&
139
         &          zat,input%xa,input%thetad,input%epsdisp,input%epsforce,tote,&
140
         &          xold,y,h,tau0, lconv)
141

142 143
    IF (lconv) THEN
       WRITE (6,'(a)') "Des woars!"
144
       CALL juDFT_end(" GEO Des woars ", 0) ! The 0 is temporarily. Should be mpi%irank.
145 146 147
    ELSE
       na = 0
       DO itype=1,atoms%ntype
148
          tau0_i(:,itype)=MATMUL(cell%bmat,tau0(:,itype))/tpi_const
149 150 151 152 153 154 155 156
          DO ieq = 1,atoms%neq(itype)
             na = na + 1
             jop = sym%invtab(atoms%ngopr(na))
             IF (oneD%odi%d1) jop = oneD%ods%ngopr(na)
             DO i = 1,3
                atoms_new%taual(i,na) = 0.0
                DO j = 1,3
                   IF (.NOT.oneD%odi%d1) THEN
157
                      atoms_new%taual(i,na) = atoms_new%taual(i,na) + sym%mrot(i,j,jop) * tau0_i(j,itype)
158
                   ELSE
159
                      atoms_new%taual(i,na) = atoms_new%taual(i,na) + oneD%ods%mrot(i,j,jop) * tau0_i(j,itype)
160 161 162
                   END IF
                ENDDO
                IF (oneD%odi%d1) THEN
163
                   atoms_new%taual(i,na) = atoms_new%taual(i,na) + oneD%ods%tau(i,jop)/cell%amat(3,3)
164 165 166 167 168 169 170
                ELSE
                   atoms_new%taual(i,na) = atoms_new%taual(i,na) + sym%tau(i,jop)
                END IF
             ENDDO
          ENDDO
       ENDDO

171 172
       input%l_f = .FALSE.

173
       IF(.NOT.input%l_inpXML) THEN
174
          CALL judft_error("Relaxation no longer supported with old inp-file") 
175
       ELSE
176
          kpts_temp%ntet = 1
177 178
          kpts_temp%numSpecialPoints = 1
          ALLOCATE(kpts_temp%specialPoints(3,kpts_temp%numSpecialPoints))
179 180 181
          ALLOCATE(noel_temp(1),atomTypeSpecies(1),speciesRepAtomType(1))
          ALLOCATE(xmlElectronStates(1,1),xmlPrintCoreStates(1,1))
          ALLOCATE(xmlCoreOccs(1,1,1))
182
          CALL initWannierDefaults(wann_temp)
183
          CALL r_inpXML(atoms_temp,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,&
184
                        banddos_temp,dimension_temp,forcetheo,cell_temp,sym_temp,xcpot_temp,noco_temp,&
185
                        oneD_temp,hybrid_temp,kpts_temp,enpara_temp,coreSpecInput_temp,wann_temp,noel_temp,&
186
                        namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,dtild_temp,xmlElectronStates,&
187
                        xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,l_kpts_temp)
188 189
          numSpecies = SIZE(speciesRepAtomType)
          filename = 'inp_new.xml'
190 191
          input_temp = input
          !input_temp%l_f = input%l_f
192
          input_temp%gw_neigd = dimension_temp%neigd
Daniel Wortmann's avatar
Daniel Wortmann committed
193
          div(:) = MIN(kpts_temp%nkpt3(:),1)
194
          stars_temp%gmax = stars_temp%gmaxInit
195
          CALL w_inpXML(atoms_new,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,forcetheo,&
196
                        banddos_temp,cell_temp,sym_temp,xcpot_temp,noco_temp,oneD_temp,hybrid_temp,&
Daniel Wortmann's avatar
Daniel Wortmann committed
197
                        kpts_temp,kpts_temp%nkpt3,kpts_temp%l_gamma,noel_temp,namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,&
198
                        dtild_temp,input_temp%comment,xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs,&
199
                        atomTypeSpecies,speciesRepAtomType,.FALSE.,filename,.TRUE.,numSpecies,enpara_temp)
200
          DEALLOCATE(atomTypeSpecies,speciesRepAtomType)
201 202
          DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs)
       END IF
203

Daniel Wortmann's avatar
Daniel Wortmann committed
204
    ENDIF
205 206 207
    RETURN
  END SUBROUTINE geo
END MODULE m_geo