geo.f90 8.66 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 10 11 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 41
MODULE m_geo
  USE m_juDFT
CONTAINS
  SUBROUTINE geo(&
       &               atoms,sym,cell,oneD,input_in,tote,&
       &               forcetot)

    !    *********************************************************************
    !    * 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
42 43
    USE m_rinpXML
    USE m_winpXML
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    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
    TYPE(t_input),INTENT(IN):: input_in
    ! ..
    ! ..  Scalar Arguments ..
    REAL,    INTENT (IN) :: tote
    ! ..
    ! ..  Array Arguments ..
    REAL,    INTENT (INOUT) :: forcetot(3,atoms%ntypd)
    ! ..
    ! ..  Local Scalars ..
    INTEGER i,j,na ,istep0,istep,itype,jop,ieq
    LOGICAL lconv       
    TYPE(t_atoms)  :: atoms_new
    ! ..
    ! ..  Local Arrays ..
    REAL xold(3*atoms%ntypd),y(3*atoms%ntypd),h(3*atoms%ntypd,3*atoms%ntypd),zat(atoms%ntypd)
    REAL tau0(3,atoms%ntypd),tau0_i(3,atoms%ntypd) 

    TYPE(t_input):: input

68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
    ! temporary variables for XML IO
    TYPE(t_input)                 :: input_temp
    TYPE(t_dimension)             :: dimension_temp
    TYPE(t_atoms)                 :: atoms_temp
    TYPE(t_sphhar)                :: sphhar_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
    LOGICAL                       :: l_opti_temp
    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

103
    input=input_in
104 105
    atoms_new=atoms

106 107 108 109 110
    na = 1
    DO i = 1,atoms_new%ntype
       IF (input%film) atoms_new%taual(3,na) = atoms_new%taual(3,na)/cell%amat(3,3)
       tau0_i(:,i) = atoms_new%taual(:,na)
       tau0(:,i)=MATMUL(cell%amat,tau0_i(:,i))
111
       tau0(:,i)=atoms%pos(:,na)
112 113
       na = na + atoms_new%neq(i)
    END DO
114

115

116
    CALL bfgs0(atoms%ntype, istep0,xold,y,h)
117 118 119

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

    istep = 1
134
    CALL bfgs(atoms%ntype,istep,istep0,forcetot,&
135
         &          zat,input%xa,input%thetad,input%epsdisp,input%epsforce,tote,&
136
         &          xold,y,h,tau0, lconv)
137 138
    IF (lconv) THEN
       WRITE (6,'(a)') "Des woars!"
139
       CALL juDFT_end(" GEO Des woars ", 1) ! The 1 is temporarily. Should be mpi%irank.
140 141 142 143 144 145 146 147 148 149 150 151
    ELSE
       na = 0
       DO itype=1,atoms%ntype
          tau0_i(:,itype)=MATMUL(cell%bmat,tau0(:,itype))
          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
152
                      atoms_new%taual(i,na) = atoms_new%taual(i,na) + sym%mrot(i,j,jop) * tau0_i(j,itype)
153
                   ELSE
154
                      atoms_new%taual(i,na) = atoms_new%taual(i,na) + oneD%ods%mrot(i,j,jop) * tau0_i(j,itype)
155 156 157
                   END IF
                ENDDO
                IF (oneD%odi%d1) THEN
158
                   atoms_new%taual(i,na) = atoms_new%taual(i,na) + oneD%ods%tau(i,jop)/cell%amat(3,3)
159 160 161 162 163 164 165
                ELSE
                   atoms_new%taual(i,na) = atoms_new%taual(i,na) + sym%tau(i,jop)
                END IF
             ENDDO
          ENDDO
       ENDDO

166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
!       CALL judft_error("Writing on new input file not implemented in geo")
       input%l_f = .false.
!       CALL rw_inp('W',atoms_new,obsolete,vacuum,input,stars,sliceplot,banddos,&
!                   cell,sym,xcpot,noco,jij,oneD,hybrid,kpts,&
!                   noel,namex,relcor,a1,a2,a3,scale,dtild,name)
       IF(input%l_inpXML) THEN
          ALLOCATE(noel_temp(1),atomTypeSpecies(1),speciesRepAtomType(1))
          ALLOCATE(xmlElectronStates(1,1),xmlPrintCoreStates(1,1))
          ALLOCATE(xmlCoreOccs(1,1,1))
          CALL r_inpXML(&
                        atoms_temp,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,&
                        banddos_temp,dimension_temp,cell_temp,sym_temp,xcpot_temp,noco_temp,Jij_temp,&
                        oneD_temp,hybrid_temp,kpts_temp,enpara_temp,sphhar_temp,l_opti_temp,noel_temp,&
                        namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,scale_temp,dtild_temp,xmlElectronStates,&
                        xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType)
          numSpecies = SIZE(speciesRepAtomType)
          filename = 'inp_new.xml'
          input_temp%l_f = input%l_f
          div(:) = MIN(kpts_temp%nmop(:),1)
          CALL w_inpXML(&
186
                        atoms_new,obsolete_temp,vacuum_temp,input_temp,stars_temp,sliceplot_temp,&
187 188 189
                        banddos_temp,cell_temp,sym_temp,xcpot_temp,noco_temp,jij_temp,oneD_temp,hybrid_temp,&
                        kpts_temp,kpts_temp%nmop,kpts_temp%l_gamma,noel_temp,namex_temp,relcor_temp,a1_temp,a2_temp,a3_temp,&
                        scale_temp,dtild_temp,input_temp%comment,xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs,&
190
                        atomTypeSpecies,speciesRepAtomType,.FALSE.,filename,numSpecies,enpara_temp)
191 192 193
          DEALLOCATE(noel_temp,atomTypeSpecies,speciesRepAtomType)
          DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs)
       END IF
194 195 196 197 198
    ENDIF

    RETURN
  END SUBROUTINE geo
END MODULE m_geo