geo.f90 9.83 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
    USE m_xsf_io
46 47 48 49 50
    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
51 52
    TYPE(t_vacuum),INTENT(IN) :: vacuum
    TYPE(t_input),INTENT(IN)  :: input_in
53 54 55 56 57
    ! ..
    ! ..  Scalar Arguments ..
    REAL,    INTENT (IN) :: tote
    ! ..
    ! ..  Array Arguments ..
58
    REAL,    INTENT (INOUT) :: forcetot(3,atoms%ntype)
59 60 61
    ! ..
    ! ..  Local Scalars ..
    INTEGER i,j,na ,istep0,istep,itype,jop,ieq
62
    LOGICAL lconv
63 64 65
    TYPE(t_atoms)  :: atoms_new
    ! ..
    ! ..  Local Arrays ..
66 67
    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) 
68 69 70

    TYPE(t_input):: input

71 72 73 74 75 76 77 78 79 80 81 82 83
    ! 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
84
    CLASS(t_xcpot),ALLOCATABLE    :: xcpot_temp
85 86 87 88
    TYPE(t_results)               :: results_temp
    TYPE(t_kpts)                  :: kpts_temp
    TYPE(t_hybrid)                :: hybrid_temp
    TYPE(t_oneD)                  :: oneD_temp
89
    TYPE(t_coreSpecInput)         :: coreSpecInput_temp
90
    TYPE(t_wann)                  :: wann_temp
91
    LOGICAL                       :: l_kpts_temp
92 93 94 95 96 97 98 99 100 101 102 103 104
    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
105
    REAL                          :: forceAllAtoms(3,atoms%nat)
106
    CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
107
    TYPE(t_field)                 :: field
108
    input=input_in
109 110
    atoms_new=atoms

111 112 113 114 115
    istep0 = 0
    xold = 0.0
    y = 0.0
    h= 0.0

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

123
    CALL bfgs0(atoms%ntype, istep0,xold,y,h)
124 125 126

    DO itype=1,atoms%ntype
       IF (atoms%l_geo(itype)) THEN
127
          WRITE (6,'(6f10.5)') (tau0(j,itype),j=1,3), (forcetot(i,itype),i=1,3)
128 129 130 131
          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),&
132
               &     (forcetot(i,itype),i=1,3),' atoms%relax: ', (atoms%relax(i,itype),i=1,3)
133 134 135 136 137 138 139 140
       ELSE
          DO i = 1,3
             forcetot(i,itype)=0.0
          ENDDO
       ENDIF
    ENDDO

    istep = 1
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
    CALL bfgs(atoms%ntype,istep,istep0,forcetot,zat,input%xa,input%thetad,input%epsdisp,&
              input%epsforce,tote,xold,y,h,tau0,lconv)

    !write out struct_force.xsf file
    forceAllAtoms = 0.0
    na = 0
    DO itype=1,atoms%ntype
       forcetot(:,itype)=MATMUL(cell%bmat,forcetot(:,itype))/tpi_const ! to inner coordinates
       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
             DO j = 1,3
                IF (.NOT.oneD%odi%d1) THEN
                   forceAllAtoms(i,na) = forceAllAtoms(i,na) + sym%mrot(i,j,jop) * forcetot(j,itype)
                ELSE
                   forceAllAtoms(i,na) = forceAllAtoms(i,na) + oneD%ods%mrot(i,j,jop) * forcetot(j,itype)
                END IF
             END DO
          END DO
          forceAllAtoms(:,na) = MATMUL(cell%amat,forceAllAtoms(:,na)) ! to external coordinates
       END DO
    END DO
    OPEN (55,file="struct_force.xsf",status='replace')
    CALL xsf_WRITE_atoms(55,atoms,input%film,.false.,cell%amat,forceAllAtoms)
    CLOSE (55)
168

169 170
    IF (lconv) THEN
       WRITE (6,'(a)') "Des woars!"
171
       CALL juDFT_end(" GEO Des woars ", 0) ! The 0 is temporarily. Should be mpi%irank.
172 173 174
    ELSE
       na = 0
       DO itype=1,atoms%ntype
175
          tau0_i(:,itype)=MATMUL(cell%bmat,tau0(:,itype))/tpi_const
176 177 178 179 180 181 182 183
          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
184
                      atoms_new%taual(i,na) = atoms_new%taual(i,na) + sym%mrot(i,j,jop) * tau0_i(j,itype)
185
                   ELSE
186
                      atoms_new%taual(i,na) = atoms_new%taual(i,na) + oneD%ods%mrot(i,j,jop) * tau0_i(j,itype)
187 188 189
                   END IF
                ENDDO
                IF (oneD%odi%d1) THEN
190
                   atoms_new%taual(i,na) = atoms_new%taual(i,na) + oneD%ods%tau(i,jop)/cell%amat(3,3)
191 192 193 194 195 196 197
                ELSE
                   atoms_new%taual(i,na) = atoms_new%taual(i,na) + sym%tau(i,jop)
                END IF
             ENDDO
          ENDDO
       ENDDO

198 199
       input%l_f = .FALSE.

200
       IF(.NOT.input%l_inpXML) THEN
201
          CALL judft_error("Relaxation no longer supported with old inp-file") 
202
       ELSE
203
          kpts_temp%ntet = 1
204 205
          kpts_temp%numSpecialPoints = 1
          ALLOCATE(kpts_temp%specialPoints(3,kpts_temp%numSpecialPoints))
206 207 208
          ALLOCATE(noel_temp(1),atomTypeSpecies(1),speciesRepAtomType(1))
          ALLOCATE(xmlElectronStates(1,1),xmlPrintCoreStates(1,1))
          ALLOCATE(xmlCoreOccs(1,1,1))
209
          CALL initWannierDefaults(wann_temp)
210
          CALL r_inpXML(atoms_temp,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,&
211
                        banddos_temp,dimension_temp,forcetheo,field,cell_temp,sym_temp,xcpot_temp,noco_temp,&
212
                        oneD_temp,hybrid_temp,kpts_temp,enpara_temp,coreSpecInput_temp,wann_temp,noel_temp,&
213
                        namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,dtild_temp,xmlElectronStates,&
214
                        xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,l_kpts_temp)
215 216
          numSpecies = SIZE(speciesRepAtomType)
          filename = 'inp_new.xml'
217 218
          input_temp = input
          !input_temp%l_f = input%l_f
219
          input_temp%gw_neigd = dimension_temp%neigd
Daniel Wortmann's avatar
Daniel Wortmann committed
220
          div(:) = MIN(kpts_temp%nkpt3(:),1)
221
          stars_temp%gmax = stars_temp%gmaxInit
222
          CALL w_inpXML(atoms_new,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,forcetheo,&
223
                        banddos_temp,cell_temp,sym_temp,xcpot_temp,noco_temp,oneD_temp,hybrid_temp,&
Daniel Wortmann's avatar
Daniel Wortmann committed
224
                        kpts_temp,kpts_temp%nkpt3,kpts_temp%l_gamma,noel_temp,namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,&
225
                        dtild_temp,input_temp%comment,xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs,&
226
                        atomTypeSpecies,speciesRepAtomType,.FALSE.,filename,.TRUE.,numSpecies,enpara_temp)
227
          DEALLOCATE(atomTypeSpecies,speciesRepAtomType)
228 229
          DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs)
       END IF
230

Daniel Wortmann's avatar
Daniel Wortmann committed
231
    ENDIF
232 233 234
    RETURN
  END SUBROUTINE geo
END MODULE m_geo