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 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

    !    *********************************************************************
    !    * calculates the NEW atomic positions after the results%force calculation   *
    !    * SUBROUTINE is based on a BFGS method implemented by jij%M. Weinert    *
    !    *                                [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
    !
    !    *********************************************************************
    USE m_rwinp
    USE m_bfgs
    USE m_bfgs0
    USE m_types
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 83 84 85 86 87 88
    ! 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
    TYPE(t_xcpot)                 :: xcpot_temp
    TYPE(t_results)               :: results_temp
    TYPE(t_jij)                   :: jij_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 105
    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

106
    input=input_in
107 108
    atoms_new=atoms

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

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

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

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

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

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

172 173
       input%l_f = .FALSE.

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

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