1! 2! Copyright (C) 2003-2006 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!-------------------------------------------------------------------------- 9MODULE fcp_opt_routines 10 !--------------------------------------------------------------------------- 11 ! 12 ! ... This module contains all subroutines and functions needed for 13 ! ... the optimisation of the reaction path (NEB and SMD calculations) 14 ! 15 ! ... Written by Carlo Sbraccia ( 2003-2006 ) 16 ! 17 ! ... MDIIS algorithm is implemented by Satomichi Nishihara ( 2016 ) 18 ! 19 USE kinds, ONLY : DP 20 USE constants, ONLY : eps8, eps16, e2, rytoev, fpi 21 USE path_variables, ONLY : ds, pos, grad 22 USE io_global, ONLY : meta_ionode, meta_ionode_id 23 USE mp, ONLY : mp_bcast 24 USE mp_world, ONLY : world_comm 25 USE fcp_variables, ONLY : fcp_mu, fcp_relax, fcp_relax_step, & 26 fcp_mdiis_size, fcp_mdiis_step, & 27 fcp_tot_charge_first, fcp_tot_charge_last 28 USE mdiis, ONLY : mdiis_type, allocate_mdiis, deallocate_mdiis, update_by_mdiis 29 USE path_variables, ONLY : num_of_images 30 ! 31 IMPLICIT NONE 32 ! 33 PRIVATE 34 ! 35 REAL(DP), ALLOCATABLE :: fcp_neb_nelec(:) 36 REAL(DP), ALLOCATABLE :: fcp_neb_ef(:) 37 ! 38 ! ... variables for line-minimisation 39 REAL(DP), ALLOCATABLE :: force0(:) 40 REAL(DP), ALLOCATABLE :: nelec0(:) 41 LOGICAL, ALLOCATABLE :: firstcall(:) 42 ! 43 ! ... variables for MDIIS 44 LOGICAL :: init_mdiis 45 TYPE(mdiis_type) :: mdiist 46 ! 47 PUBLIC :: fcp_neb_nelec, fcp_neb_ef, & 48 fcp_opt_allocation, fcp_opt_deallocation, fcp_opt_perform 49 ! 50CONTAINS 51 ! 52 !---------------------------------------------------------------------- 53 SUBROUTINE fcp_opt_allocation() 54 !---------------------------------------------------------------------- 55 ! 56 USE ions_base, ONLY : nat, ityp, zv 57 USE klist, ONLY : tot_charge, nelec 58 USE io_files, ONLY : prefix, tmp_dir 59 USE path_variables, ONLY : restart 60 USE ener, ONLY : ef 61 ! 62 IMPLICIT NONE 63 ! 64 REAL(DP) :: ionic_charge, nelec_, first, last 65 INTEGER :: n, i, ierr 66 CHARACTER (LEN=256) :: tmp_dir_saved 67 ! 68 CHARACTER(LEN=6), EXTERNAL :: int_to_char 69 ! 70 ALLOCATE( fcp_neb_nelec( num_of_images ) ) 71 ALLOCATE( fcp_neb_ef ( num_of_images ) ) 72 ! 73 IF ( TRIM(fcp_relax) == 'lm' ) THEN 74 ! 75 ALLOCATE( force0 ( num_of_images ) ) 76 ALLOCATE( nelec0 ( num_of_images ) ) 77 ALLOCATE( firstcall ( num_of_images ) ) 78 ! 79 force0 (:) = 0.0_DP 80 nelec0 (:) = 0.0_DP 81 firstcall (:) = .TRUE. 82 ! 83 ELSE IF ( TRIM(fcp_relax) == 'mdiis' ) THEN 84 ! 85 init_mdiis = .TRUE. 86 CALL allocate_mdiis(mdiist, fcp_mdiis_size, num_of_images, fcp_mdiis_step, 1) 87 ! 88 END IF 89 ! 90 IF ( restart ) THEN 91 ! 92 tmp_dir_saved = tmp_dir 93 ! 94 DO i = 1, num_of_images 95 ! 96 tmp_dir = TRIM( tmp_dir_saved ) // TRIM( prefix ) // "_" // & 97 TRIM( int_to_char( i ) ) // "/" 98 ! 99 CALL errore('fcp_opt_routines','XSD implementation pending',1) 100 ! 101 fcp_neb_nelec(i) = nelec 102 fcp_neb_ef (i) = ef * e2 ! factor e2: hartree -> Ry. 103 ! 104 END DO 105 ! 106 tmp_dir = tmp_dir_saved 107 ! 108 ELSE 109 ! 110 ionic_charge = SUM(zv(ityp(1:nat))) 111 nelec_ = ionic_charge - tot_charge 112 ! 113 n = num_of_images 114 first = fcp_tot_charge_first 115 last = fcp_tot_charge_last 116 DO i = 1, n 117 fcp_neb_nelec(i) = nelec_ - (first * (n - i) + last * (i - 1) ) / (n - 1) 118 END DO 119 ! 120 fcp_neb_ef(:) = 0.0_DP 121 ! 122 END IF 123 ! 124 END SUBROUTINE fcp_opt_allocation 125 ! 126 !---------------------------------------------------------------------- 127 SUBROUTINE fcp_opt_deallocation() 128 !---------------------------------------------------------------------- 129 ! 130 IMPLICIT NONE 131 ! 132 IF ( ALLOCATED( fcp_neb_nelec ) ) DEALLOCATE( fcp_neb_nelec ) 133 IF ( ALLOCATED( fcp_neb_ef ) ) DEALLOCATE( fcp_neb_ef ) 134 IF ( ALLOCATED( force0 ) ) DEALLOCATE( force0 ) 135 IF ( ALLOCATED( nelec0 ) ) DEALLOCATE( nelec0 ) 136 IF ( ALLOCATED( firstcall ) ) DEALLOCATE( firstcall ) 137 ! 138 IF ( init_mdiis ) THEN 139 ! 140 CALL deallocate_mdiis(mdiist) 141 ! 142 END IF 143 ! 144 END SUBROUTINE fcp_opt_deallocation 145 ! 146 !---------------------------------------------------------------------- 147 SUBROUTINE fcp_opt_perform() 148 !---------------------------------------------------------------------- 149 ! 150 IMPLICIT NONE 151 ! 152 IF ( TRIM(fcp_relax) == 'lm' ) THEN 153 ! 154 CALL fcp_line_minimisation() 155 ! 156 ELSE IF ( TRIM(fcp_relax) == 'mdiis' ) THEN 157 ! 158 CALL fcp_mdiis() 159 ! 160 END IF 161 ! 162 END SUBROUTINE fcp_opt_perform 163 ! 164 !---------------------------------------------------------------------- 165 SUBROUTINE fcp_line_minimisation() 166 !---------------------------------------------------------------------- 167 ! 168 USE ions_base, ONLY : nat, ityp, zv 169 USE cell_base, ONLY : at, alat 170 ! 171 USE path_variables, ONLY : frozen 172 ! 173 IMPLICIT NONE 174 ! 175 INTEGER :: image 176 REAL(DP) :: force, ef, nelec, n_tmp, max_tot_charge, capacitance, ionic_charge 177 ! 178 IF ( meta_ionode ) THEN 179 ! 180 DO image = 1, num_of_images 181 ! 182 IF ( frozen(image) ) CYCLE 183 ! 184 ef = fcp_neb_ef(image) 185 nelec = fcp_neb_nelec(image) 186 ! 187 force = fcp_mu - ef 188 ! 189 ! ... assumption: capacitance with vacuum gives the upper bound of 190 ! ... tot_charge difference. 191 ! 192 capacitance = (at(1,1) * at(2,2) - at(2,1) * at(1,2)) * alat**2 & 193 / (alat * at(3,3) / 2._DP ) / fpi 194 max_tot_charge = abs( capacitance * force / e2 ) 195 IF ( firstcall(image) .OR. ABS( force0(image) - force ) < 1.0D-20 ) THEN 196 firstcall(image) = .FALSE. 197 nelec0(image) = nelec 198 force0(image) = force 199 nelec = nelec + fcp_relax_step * force 200 ELSE 201 n_tmp = nelec 202 nelec = (nelec * force0(image) - nelec0(image) * force ) & 203 / ( force0(image) - force ) 204 nelec0(image) = n_tmp 205 force0(image) = force 206 END IF 207 ionic_charge = SUM(zv(ityp(1:nat))) 208 !write( *,'(/,5X,"Upper bound for tot_charge:",F12.6)') & 209 ! max_tot_charge 210 !write( *,'(5X,"Original:",F12.6,"Expected:",F12.6)') & 211 ! ionic_charge - nelec0(image), ionic_charge - nelec 212 if( nelec-nelec0(image) < -max_tot_charge ) & 213 nelec= nelec0(image) - max_tot_charge 214 if( nelec-nelec0(image) > max_tot_charge ) & 215 nelec= nelec0(image) + max_tot_charge 216 !write( *,'(5X,"Next tot_charge:",F12.6)') ionic_charge - nelec 217 ! 218 fcp_neb_nelec(image) = nelec 219 ! 220 END DO 221 ! 222 END IF 223 ! 224 CALL mp_bcast( fcp_neb_nelec, meta_ionode_id, world_comm ) 225 ! 226 RETURN 227 ! 228 END SUBROUTINE fcp_line_minimisation 229 ! 230 !---------------------------------------------------------------------- 231 SUBROUTINE fcp_mdiis() 232 !---------------------------------------------------------------------- 233 ! 234 USE path_variables, ONLY : frozen 235 ! 236 IMPLICIT NONE 237 ! 238 INTEGER :: image 239 REAL(DP) :: ef 240 REAL(DP), ALLOCATABLE :: force1(:) 241 REAL(DP), ALLOCATABLE :: nelec1(:) 242 ! 243 ALLOCATE(force1(num_of_images)) 244 ALLOCATE(nelec1(num_of_images)) 245 ! 246 IF ( meta_ionode ) THEN 247 ! 248 DO image = 1, num_of_images 249 ! 250 ef = fcp_neb_ef(image) 251 nelec1(image) = fcp_neb_nelec(image) 252 force1(image) = fcp_mu - ef 253 ! 254 END DO 255 ! 256 CALL update_by_mdiis(mdiist, nelec1, force1) 257 ! 258 DO image = 1, num_of_images 259 ! 260 IF ( frozen(image) ) CYCLE 261 ! 262 fcp_neb_nelec(image) = nelec1(image) 263 ! 264 END DO 265 ! 266 END IF 267 ! 268 CALL mp_bcast( fcp_neb_nelec, meta_ionode_id, world_comm ) 269 ! 270 DEALLOCATE(force1) 271 DEALLOCATE(nelec1) 272 ! 273 RETURN 274 ! 275 END SUBROUTINE fcp_mdiis 276 ! 277END MODULE fcp_opt_routines 278