TABLE OF CONTENTS
ABINIT/scfcge [ Functions ]
NAME
scfcge
FUNCTION
Compute the next vtrial of the SCF cycle. Uses a conjugate gradient minimization of the total energy Can move only the trial potential (if moved_atm_inside==0), or move the trial atomic positions as well (if moved_atm_inside==1).
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX dtn_pc(3,natom)=preconditioned change of atomic position, in reduced coordinates. Will be quickly transferred to f_atm(:,:,i_vrespc(1)) etotal=the actual total energy initialized= if 0, the initialization of the gstate run is not yet finished iscf =5 => SCF cycle, CG based on estimation of energy gradient =6 => SCF cycle, CG based on true minimization of the energy isecur=level of security of the computation istep= number of the step in the SCF cycle moved_atm_inside: if==1, the atoms are allowed to move. mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT natom=number of atoms nfft=(effective) number of FFT grid points (for this processor) nfftot=total number of FFT grid points nspden=number of spin-density components n_fftgr=third dimension of the array f_fftgr n_index=dimension for indices of potential/density (see i_vresid, ivrespc, i_rhor...) opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential 1 vtrial (and also f_fftgr) actually contains the trial density response= if 0, GS calculation, if 1, RF calculation, intrinsically harmonic ! rhor(cplex*nfft,nspden)=actual density ucvol=unit cell volume in bohr**3
OUTPUT
dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
SIDE EFFECTS
Input/Output:
vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
the input residual of the potential and Hellman-Feynman forces
at output, it is the new trial potential .
xred(3,natom)=(needed if moved_atm_inside==1)
reduced dimensionless atomic coordinates
at input, those that generated the input residual of the potential
and Hellman-Feynman forces, at output, these are the new ones.
f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
The input vtrial is transferred, at output, in f_fftgr(:,:,1).
The input f_fftgr(:,:,i_vresid(1)) contains the last residual.
the value of i_vresid(1) is transferred to i_vresid(2) at output.
The input f_fftgr(:,:,i_vresid(2)) contains the old residual.
the value of i_vresid(2) is transferred to i_vresid(3) at output.
The input f_fftgr(:,:,i_vresid(3)) contains the previous last residual.
For the preconditioned potential residual, the same logic as for the
the potential residual is used, with i_vrespc replacing i_vresid.
The input rhor is transferred, at output, in f_fft(:,:,i_rhor(2)).
The old density is input in f_fft(:,:,i_rhor(2)), and the value of
i_rhor(2) is transferred to i_rhor(3) before the end of the routine.
The input/output search vector is stored in f_fftgr(:,:,6)
f_atm(3,natom,n_fftgr)=different functions defined for each atom :
The input xred is transferred, at output, in f_atm(:,:,1).
The input f_atm(:,:,i_vresid(1)) contains minus the HF forces.
the value of i_vresid(1) is transferred to i_vresid(2) at output.
The input f_atm(:,:,i_vresid(2)) contains minus the old HF forces.
the value of i_vresid(2) is transferred to i_vresid(3) at output.
The input f_atm(:,:,i_vresid(3)) contains minus the previous old HF forces.
For the preconditioned change of atomic positions, the same logic as for the
the potential residual is used, with i_vrespc replacing i_vresid.
The input/output search vector is stored in f_atm(:,:,6)
i_rhor(2:3)=index of the density (past and previous past) in the array f_fftgr
i_vresid(3)=index of the residual potentials (present, past and previous
past) in the array f_fftgr; also similar index for minus Hellman-Feynman
forces in the array f_atm .
i_vrespc(3)=index of the preconditioned residual potentials
(present, past and previous past) in the array f_fftgr ;
also similar index for the preconditioned change of atomic position (dtn_pc).
TODO
This routine is much too difficult to read ! Should be rewritten ... Maybe make separate subroutines for line search and CG step ?!
PARENTS
m_ab7_mixing
CHILDREN
aprxdr,findminscf,sqnormm_v,wrtout
SOURCE
101 #if defined HAVE_CONFIG_H 102 #include "config.h" 103 #endif 104 105 #include "abi_common.h" 106 107 108 subroutine scfcge(cplex,dbl_nnsclo,dtn_pc,etotal,f_atm,& 109 & f_fftgr,initialized,iscf,isecur,istep,& 110 & i_rhor,i_vresid,i_vrespc,moved_atm_inside,mpicomm,mpi_summarize,& 111 & natom,nfft,nfftot,nspden,n_fftgr,n_index,opt_denpot,response,rhor,ucvol,vtrial,xred,errid,errmess) 112 113 use defs_basis 114 use m_profiling_abi 115 use m_errors 116 117 !This section has been created automatically by the script Abilint (TD). 118 !Do not modify the following lines by hand. 119 #undef ABI_FUNC 120 #define ABI_FUNC 'scfcge' 121 use interfaces_14_hidewrite 122 use interfaces_56_mixing, except_this_one => scfcge 123 !End of the abilint section 124 125 implicit none 126 127 !Arguments ------------------------------------ 128 !scalars 129 integer,intent(in) :: cplex,initialized,iscf,isecur,istep,moved_atm_inside,mpicomm 130 integer,intent(in) :: n_fftgr,n_index,natom,nfft,nfftot,nspden,opt_denpot,response 131 integer,intent(out) :: dbl_nnsclo, errid 132 character(len = 500), intent(out) :: errmess 133 logical, intent(in) :: mpi_summarize 134 real(dp),intent(in) :: etotal,ucvol 135 !arrays 136 integer,intent(inout) :: i_rhor(n_index),i_vresid(n_index),i_vrespc(n_index) 137 real(dp),intent(in) :: dtn_pc(3,natom),rhor(cplex*nfft,nspden) 138 real(dp),intent(inout) :: f_atm(3,natom,n_fftgr) 139 real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr) 140 real(dp),intent(inout) :: vtrial(cplex*nfft,nspden),xred(3,natom) 141 142 !Local variables------------------------------- 143 !mlinmin gives the maximum number of steps in the line minimization 144 ! after which the algorithm is restarted (with a decrease of the 145 ! adaptative trial step length). This number should not be large, 146 ! since if the potential landscape is harmonic, the number of 147 ! search steps should be small. If it is large, we are not in the 148 ! harmonic region, and the CG algorithm will not be really useful, 149 ! so one can just restart the algorithm ... 150 !scalars 151 integer,parameter :: mlinmin=5 152 integer,save :: end_linmin,iline_cge,ilinear,ilinmin,isecur_eff,nlinear 153 integer,save :: number_of_restart,status 154 integer :: choice,iatom,idir,ifft,iline_cge_input,ilinmin_input,isp 155 integer :: testcg,tmp,errid_ 156 real(dp),save :: d2edv2_old2,d_lambda_old2,dedv_old2,etotal_old 157 real(dp),save :: etotal_previous,lambda_adapt,lambda_new,lambda_old,resid_old 158 real(dp) :: d2e11,d2e12,d2e22,d2edv2_new,d2edv2_old 159 real(dp) :: d2edv2_predict,d_lambda,de1,de2,dedv_mix 160 real(dp) :: dedv_new,dedv_old,dedv_predict,determ,etotal_input 161 real(dp) :: etotal_predict,gamma,lambda_input,lambda_predict2 162 real(dp) :: lambda_predict=1.0_dp,ratio,reduction 163 real(dp) :: resid_input,temp 164 character(len=500) :: message 165 !arrays 166 real(dp) :: resid_new(1) 167 real(dp), allocatable :: tmp_fft1(:,:) 168 169 ! ************************************************************************* 170 171 errid = AB7_NO_ERROR 172 dbl_nnsclo = 0 173 174 !reduction gives the level of reduction of the error in 175 !the line minimization to be reached for the minimization to be 176 !considered successfull 177 reduction=0.1_dp 178 179 !nlinear increases with the number of times the 2D minimization succeded 180 !to reach the true minimum directly. It is a measure of the 181 !degree of parabolicity of the problem, and is used to 182 !skip some steps by performing extrapolation. 183 if(istep==1)then 184 185 ! Skipping some steps is sometimes unsecure, so it is possible 186 ! to make nlinear start at a negative value - if isecur is positive 187 isecur_eff=isecur 188 nlinear=min(-isecur_eff,0) 189 ilinear=0 190 191 ! Response function calculation are intrinsically harmonic, so one 192 ! can shift isecur (by -2), and start with a positive nlinear 193 if(response==1)then 194 isecur_eff=isecur-2 195 nlinear=-isecur_eff 196 ilinear=nlinear 197 end if 198 199 iline_cge=0 200 ilinmin=0 201 end if 202 203 !Compute actual residual resid_new (residual of f_fftgr(:,:,i_vrespc(1)) 204 call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr) 205 206 !Save input residual and ilinmin for final printing 207 resid_input=resid_new(1) 208 etotal_input=etotal 209 ilinmin_input=ilinmin 210 iline_cge_input=iline_cge 211 !Transfer dtn_pc in f_atm 212 if(moved_atm_inside==1)then 213 f_atm(:,:,i_vrespc(1))=dtn_pc(:,:) 214 end if 215 216 !======================================================================= 217 !Now the routine is decomposed in three mutually exclusive parts : 218 !if(istep==1)then initialize the algorithm 219 !else if(ilinmin>0)then perform the line minimisation 220 !else if(ilinmin==0)then determine the new search direction (CG step) 221 !======================================================================= 222 223 224 !-------------------------------------- 225 !Here initialize the algorithm 226 if(istep==1)then 227 228 ! At the beginning of each gstate run, lambda_adapt is forced to have the 229 ! same value, that is 1.0_dp. In the other cases when istep=1 (at different 230 ! broyden steps, for example), the previously obtained 231 ! adaptive value is kept. 232 if(initialized==0)lambda_adapt=1.0_dp 233 lambda_old=0.0_dp 234 lambda_input=0.0_dp 235 number_of_restart=0 236 lambda_new=lambda_adapt 237 238 f_fftgr(:,:,1)=vtrial(:,:) 239 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 240 241 ! This copy must be written in F77, because of stack problems on the DECs 242 do isp=1,nspden 243 do ifft=1,cplex*nfft 244 f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1)) 245 end do 246 end do 247 vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6) 248 if(moved_atm_inside==1)then 249 f_atm(:,:,1)=xred(:,:) 250 f_atm(:,:,i_rhor(2))=xred(:,:) 251 ! There shouldn t be problems with the stack size for this small array. 252 f_atm(:,:,6)=f_atm(:,:,i_vrespc(1)) 253 xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6) 254 end if 255 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 256 tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 257 ilinmin=1 258 resid_old=resid_new(1) 259 etotal_old=etotal 260 261 status=0 262 263 ! -------------------------------------- 264 265 ! Here performs the line minimisation 266 else if(ilinmin>0)then 267 268 lambda_input=lambda_new 269 270 ! The choice with the Brent algorithm has been abandoned in version 1.6.m 271 272 ! Compute the approximate energy derivatives dedv_new and dedv_old, 273 ! from vresid and vresid_old 274 choice=2 275 call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,& 276 & f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,& 277 & natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred) 278 d_lambda=lambda_new-lambda_old 279 dedv_old=dedv_old/d_lambda 280 dedv_new=dedv_new/d_lambda 281 282 ! DEBUG 283 ! write(std_out,'(a,4es12.4,i3)' )' scfcge:lold,lnew,dold,dnew,status', & 284 ! & lambda_old,lambda_new,dedv_old,dedv_new,status 285 ! ENDDEBUG 286 287 if(status==0 .or. status==3)then 288 ! 289 ! Then, compute a predicted point along the line 290 ! The value of choice determines the minimization algorithm 291 ! choice=1 uses the two values of the derivative of the energy 292 ! choice=2 uses the two values of the energy, and and estimate of the 293 ! second derivative at the mid-point. 294 295 choice=1 296 if(iscf==6)choice=2 297 call findminscf(choice,dedv_new,dedv_old,dedv_predict,& 298 & d2edv2_new,d2edv2_old,d2edv2_predict,& 299 & etotal,etotal_old,etotal_predict,& 300 & lambda_new,lambda_old,lambda_predict,errid_,message) 301 if (errid_ /= AB7_NO_ERROR) then 302 call wrtout(std_out,message,'COLL') 303 end if 304 305 ! Suppress the next line for debugging (there is another such line) 306 status=0 307 308 ! DEBUG 309 ! Keep this debugging feature : it gives access to the investigation of lines 310 ! in a different approach 311 ! if(response==1 .and. istep>8)then 312 ! lambda_predict=1.2d-2 313 ! if(istep>=15)lambda_predict=lambda_predict-0.002 314 ! if(istep>=14)stop 315 ! status=3 316 ! end if 317 ! ENDDEBUG 318 319 else 320 if(status/=-1)then 321 status=-1 322 lambda_predict=-2.5_dp 323 else 324 lambda_predict=lambda_predict+0.1_dp 325 end if 326 end if 327 328 ! If the predicted point is very close to the most recent 329 ! computed point, while this is the first trial on this line, 330 ! then we are in the linear regime : 331 ! nlinear is increased by one unit. For the time being, do this even when 332 ! moved_atm_inside==1 (the code still works when it is done, but it 333 ! seems to be a bit unstable). The maximal value of nlinear is 1, except 334 ! when isecur_eff is a negative number, less than -1. 335 if( abs(lambda_predict-lambda_new)/& 336 & (abs(lambda_predict)+abs(lambda_new)) < 0.01 .and. ilinmin==1 ) then 337 ! if(moved_atm_inside==0 .and. nlinear<max(1,-isecur_eff) )nlinear=nlinear+1 338 if(nlinear<max(1,-isecur_eff))nlinear=nlinear+1 339 ilinear=nlinear 340 end if 341 342 ! If the predicted point is close to the most recent computed point, 343 ! or the previous one, set on the flag of end of line minization 344 end_linmin=0 345 if(abs(lambda_new-lambda_predict)*2.0_dp& 346 & /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1 347 if(abs(lambda_old-lambda_predict)*2.0_dp& 348 & /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1 349 350 if(status/=0)end_linmin=0 351 352 ! Save the closest old lambda, if needed, 353 ! also examine the reduction of the interval, and eventual stop 354 ! the present line minimisation, because of convergence (end_linmin=1) 355 ! Also treat the case in which the predicted value of lambda is negative, 356 ! or definitely too small in which case the algorithm has to be restarted 357 ! (not a very good solution, though ...) 358 ! Finally also treat the case where insufficiently converged 359 ! density at lambda=0.0_dp happens, which screws up the line minimisation. 360 361 ! Here restart the algorithm with the best vtrial. 362 ! Also make reduction in lambda_adapt 363 ! DEBUG 364 ! write(std_out,*)' scfcge : status=',status 365 ! ENDDEBUG 366 if( end_linmin==0 .and. status==0 .and. & 367 & ( (lambda_predict<0.005_dp*lambda_adapt .and. iscf==5) .or. & 368 & (abs(lambda_predict)<0.005_dp*lambda_adapt .and. iscf==6).or. & 369 & ilinmin==mlinmin ) )then 370 if(number_of_restart>12)then 371 errid = AB7_ERROR_MIXING_CONVERGENCE 372 write(errmess,'(a,a,i0,a,a,a,a,a)')& 373 & 'Potential-based CG line minimization not',' converged after ',number_of_restart,' restarts. ',ch10,& 374 & 'Action : read the eventual warnings about lack of convergence.',ch10,& 375 & 'Some might be relevant. Otherwise, raise nband. Returning' 376 MSG_WARNING(errmess) 377 return 378 end if 379 ! Make reduction in lambda_adapt (kind of steepest descent...) 380 write(message,'(a,a,a)')& 381 & 'Potential-based CG line minimization has trouble to converge.',ch10,& 382 & 'The algorithm is restarted with more secure parameters.' 383 MSG_WARNING(message) 384 number_of_restart=number_of_restart+1 385 ! At the second restart, double the number of non-self consistent loops. 386 if(number_of_restart>=2)dbl_nnsclo=1 387 lambda_adapt=lambda_adapt*0.7_dp 388 lambda_new=lambda_adapt 389 ! If the last energy is better than the old one, transfer the data. 390 ! Otherwise, no transfer must occur (very simple to code...) 391 if(etotal<etotal_old .or. abs(lambda_old)<1.0d-8)then 392 f_fftgr(:,:,1)=vtrial(:,:) 393 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 394 do isp=1,nspden 395 do ifft=1,cplex*nfft 396 f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1)) 397 end do 398 end do 399 if(moved_atm_inside==1)then 400 f_atm(:,:,1)=xred(:,:) 401 f_atm(:,:,i_rhor(2))=xred(:,:) 402 f_atm(:,:,6)=f_atm(:,:,i_vrespc(1)) 403 end if 404 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 405 tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 406 resid_old=resid_new(1) 407 etotal_old=etotal 408 end if 409 lambda_old=0.0_dp 410 ilinmin=1 411 ! Putting the flag to -1 avoids the usual actions taken with end_linmin=1 412 end_linmin=-1 413 ! Also put ilinear and nlinear to 0 414 ilinear=0 415 nlinear=0 416 417 ! Here lambda_new is the closest to lambda_predict, 418 ! or lambda_old is still 0.0_dp, while the energy shows that the minimum 419 ! is away from 0.0_dp (insufficiently converged density at lambda=0.0_dp). 420 else if( abs(lambda_new-lambda_predict)<abs(lambda_old-lambda_predict) & 421 & .or. & 422 & ( abs(lambda_old)<1.0d-6 .and. & 423 & ilinmin>1 .and. & 424 & etotal>etotal_previous ) & 425 & )then 426 f_fftgr(:,:,1)=vtrial(:,:) 427 tmp=i_rhor(3) ; i_rhor(3)=i_rhor(2) ; i_rhor(2)=tmp 428 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 429 tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) 430 i_vrespc(2)=i_vrespc(1); i_vrespc(1)=tmp; 431 tmp=i_vresid(3); i_vresid(3)=i_vresid(2) 432 i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 433 if(moved_atm_inside==1)then 434 f_atm(:,:,1)=xred(:,:) 435 f_atm(:,:,i_rhor(2))=xred(:,:) 436 end if 437 d_lambda_old2=lambda_old-lambda_new 438 lambda_old=lambda_new 439 etotal_old=etotal 440 resid_old=resid_new(1) 441 d2edv2_old2=d2edv2_new 442 dedv_old=dedv_new 443 dedv_old2=dedv_new 444 ! if(abs(lambda_new-lambda_predict)*2.0_dp& 445 ! & /abs(lambda_new+lambda_predict) <reduction) end_linmin=1 446 447 ! Here lambda_old is the closest to lambda_predict (except for avoiding 448 ! lambda_old==0.0_dp) 449 else 450 tmp=i_vresid(3) ; i_vresid(3)=i_vresid(1) ; i_vresid(1)=tmp 451 f_fftgr(:,:,i_rhor(3))=rhor(:,:) 452 if(moved_atm_inside==1) f_atm(:,:,i_rhor(3))=xred(:,:) 453 tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(1) ; i_vrespc(1)=tmp 454 d_lambda_old2=lambda_new-lambda_old 455 etotal_previous=etotal 456 d2edv2_old2=d2edv2_old 457 dedv_old2=dedv_old 458 ! if(abs(lambda_old-lambda_predict)*2.0_dp& 459 ! & /abs(lambda_old+lambda_predict) <reduction) end_linmin=1 460 end if 461 462 ! If the interval has not yet been sufficiently reduced, 463 ! continue the search 464 if(end_linmin==0)then 465 lambda_new=lambda_predict 466 467 ! DEBUG 468 ! write(std_out,'(a,2es16.6)' )& 469 ! & ' scfcge : continue search, lambda_old,lambda_new=',lambda_old,lambda_new 470 ! write(std_out,'(a,2es16.6)' )& 471 ! & ' scfcge : f_fftgr(3:4,1,1)=',f_fftgr(3:4,1,1) 472 ! write(std_out,'(a,2es16.6)' )& 473 ! & ' scfcge : f_fftgr(3:4,1,6)=',f_fftgr(3:4,1,6) 474 ! ENDDEBUG 475 476 vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6) 477 if(moved_atm_inside==1)then 478 xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6) 479 end if 480 481 ilinmin=ilinmin+1 482 ! 483 ! Here generates a starting point for next line search 484 else 485 iline_cge=iline_cge+1 486 if(end_linmin==1)ilinmin=0 487 lambda_old=0.0_dp 488 489 ! In order to generate the new step, take into account previous 490 ! optimal lambdas (including those of previous ion moves), 491 ! and the selected new one, if it is positive. 492 ! However, wait iline_cge>1 to select new ones. 493 ! lambda_adapt has been initialized at 1.0_dp 494 if(iline_cge>1 .and. lambda_new>0.0_dp )then 495 ! Actually compute a geometric mean 496 lambda_adapt= ( lambda_adapt**(dble(iline_cge-1)) * abs(lambda_new)) & 497 & **(1.0_dp/dble(iline_cge)) 498 ! In order to recover the previous algorithm, it is enough 499 ! to decomment the next line 500 ! lambda_adapt=1.0_dp 501 end if 502 lambda_new=lambda_adapt 503 504 vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2)) 505 if(moved_atm_inside==1)then 506 xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2)) 507 end if 508 509 ! End choice between continue line minim and determine new direction 510 end if 511 512 ! 513 ! ------------------------------- 514 515 ! Here perform the CG step 516 517 else if(ilinmin==0)then 518 519 ! Compute the approximate energy derivatives dedv_mix,dedv_new,dedv_old 520 choice=3 521 call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,& 522 & f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,& 523 & natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred) 524 525 dedv_mix=dedv_mix/lambda_new 526 dedv_new=dedv_new/lambda_new 527 dedv_old=dedv_old/lambda_new 528 529 ! DEBUG 530 ! write(message, '(a,3es12.4)' )' scfcge: lambda_adapt',& 531 ! & lambda_adapt 532 ! call wrtout(std_out,message,'COLL') 533 534 ! write(message, '(a,3es12.4)' )' scfcge: dedv_old,dedv_new,dedv_mix',& 535 ! & dedv_old,dedv_new,dedv_mix 536 ! call wrtout(std_out,message,'COLL') 537 ! ENDDEBUG 538 539 ! Then, compute a predicted point, either along the line, 540 ! or in a 2D plane 541 testcg=1 542 if(testcg==0)then 543 ! This part corresponds to steepest descent, 544 ! in which the line minimisation can be done 545 ! using different algorithms, varying with the value of choice 546 choice=1 547 if(iscf==6)choice=2 548 call findminscf(choice,dedv_new,dedv_old,dedv_predict,& 549 & d2edv2_new,d2edv2_old,d2edv2_predict,& 550 & etotal,etotal_old,etotal_predict,& 551 & lambda_new,lambda_old,lambda_predict,errid_,message) 552 if (errid_ /= AB7_NO_ERROR) then 553 call wrtout(std_out,message,'COLL') 554 end if 555 lambda_predict2=0.0_dp 556 ! Suppress the next line for debugging (there is another such line) 557 status=0 558 else 559 ! This part corresponds to conjugate gradient 560 ! A 2D minimisation is performed 561 ! oldest direction is labelled 2 562 ! newest direction is labelled 1 563 de1=dedv_old ; de2=dedv_old2 564 d2e11=(dedv_new-dedv_old)/lambda_new 565 d2e22=d2edv2_old2 566 d2e12=(dedv_mix-dedv_old)/d_lambda_old2 567 ! The system to be solved is 568 ! 0 = de1 + lambda1 d2e11 + lambda2 d2d12 569 ! 0 = de2 + lambda1 d2e12 + lambda2 d2d22 570 determ=d2e11*d2e22-d2e12*d2e12 571 lambda_predict=-(de1*d2e22-de2*d2e12)/determ 572 lambda_predict2=(de1*d2e12-de2*d2e11)/determ 573 d2edv2_new=d2e11 ; d2edv2_old=d2e11 574 end if 575 576 ! DEBUG 577 ! write(message, '(a,5es11.3)' )' scfcge: de1,de2,d2e11,d2e22,d2e12',& 578 ! & de1,de2,d2e11,d2e22,d2e12 579 ! call wrtout(std_out,message,'COLL') 580 ! write(std_out,'(a,2es12.4)' )' scfcge: la_predict,la_predict2',& 581 ! & lambda_predict,lambda_predict2 582 ! ----- 583 ! write(std_out,*)'residues ', 584 ! !$ de1+lambda_predict*d2e11+lambda_predict2*d2e12, 585 ! !$ de2+lambda_predict*d2e12+lambda_predict2*d2e22 586 ! if(.true.)stop 587 ! ENDDEBUG 588 ! 589 590 ! Determine the region of the 2D search space 591 ! in which the predicted point is located, 592 ! or use linear indicator to decide interpolation 593 ! and advance to next 2D search. 594 end_linmin=0 595 write(message, '(a,2i3)' )' nlinear, ilinear',nlinear,ilinear 596 call wrtout(std_out,message,'COLL') 597 if(lambda_predict<0.0_dp)then 598 ! Something is going wrong. Just take a reasonable step 599 ! along the steepest descent direction (Region III). 600 ! Actually, Region I and region III are treated in the same way later. 601 ! In effect, this corresponds to restart the algorithm 602 end_linmin=3 603 ! Also put ilinear and nlinear to 0 604 ilinear=0 605 nlinear=0 606 ! Decrease the adaptive step to predict next direction 607 lambda_adapt=lambda_adapt*0.7_dp 608 else if(ilinear>=1) then 609 ! Region IV : will do an interpolation 610 end_linmin=4 611 ilinear=ilinear-1 612 else if(abs(lambda_predict2)>reduction .or.& 613 & lambda_predict<0.5_dp .or.& 614 & lambda_predict>2.5_dp .or.& 615 & lambda_predict-abs(lambda_predict2)/reduction <0.0_dp ) then 616 ! Region II : lambda_predict is not too good, and not too bad. 617 end_linmin=2 618 else if (abs(1.0_dp-lambda_predict)<reduction)then 619 ! Region I, the out-of-line point is OK. 620 end_linmin=1 621 else 622 ! If everything fails, then region II. 623 end_linmin=2 624 end if 625 626 ! DEBUG 627 ! write(message, '(a,2es12.4,i2)' )& 628 ! & ' scfcge : la_predict, la_predict2, region',& 629 ! & lambda_predict,lambda_predict2,end_linmin 630 ! call wrtout(std_out,message,'COLL') 631 ! ENDDEBUG 632 633 ! Treat region I, in the same way as region III 634 if(end_linmin==1 .or. end_linmin==3)then 635 636 ! In region I, the line search is 637 ! along vtrial-vtrial_old. 638 ! The closest point is the new point 639 ! thus to be transfered in the "old" locations 640 641 do isp=1,nspden 642 do ifft=1,cplex*nfft 643 f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new 644 end do 645 end do 646 f_fftgr(:,:,1)=vtrial(:,:) 647 f_fftgr(:,:,i_rhor(2))=rhor(:,:) 648 if(moved_atm_inside==1)then 649 f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new 650 f_atm(:,:,1)=xred(:,:) 651 f_atm(:,:,i_rhor(2))=xred(:,:) 652 end if 653 tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp 654 tmp=i_vresid(3) ; i_vresid(3)=i_vresid(2) 655 i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp 656 d_lambda_old2=-lambda_new 657 lambda_old=lambda_new 658 etotal_old=etotal 659 resid_old=resid_new(1) 660 d2edv2_old=d2edv2_new 661 dedv_old=dedv_new 662 663 ! Region I or III : one is close of the 2D minimum, 664 ! or lambda_predict was negative (indicate a problem of convergence) 665 ! Compute next trial potential along the 666 ! PC residual and not along this search direction. 667 ilinmin=0 668 ! Question : isn t it here that one should prevent region I to called 669 ! itself more than 1 time ??? 670 ! Here the small difference between region I and region III 671 if(end_linmin==3)ilinmin=1 672 lambda_old=0.0_dp 673 lambda_new=lambda_adapt 674 675 vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2)) 676 if(moved_atm_inside==1)then 677 xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2)) 678 end if 679 ! The new vtrial has been generated 680 681 else 682 683 ! Here region II or IV 684 ilinmin=1 685 if (lambda_predict==0._dp) then 686 gamma=zero 687 else 688 gamma=lambda_predict2/lambda_predict 689 end if 690 ! Compute new search direction and trial potential 691 write(message,*)' compute new search direction ' 692 call wrtout(std_out,message,'COLL') 693 do isp=1,nspden 694 do ifft=1,cplex*nfft 695 f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new+ & 696 & gamma*f_fftgr(ifft,isp,6) 697 end do 698 end do 699 vtrial(:,:)=f_fftgr(:,:,1)+ lambda_predict*f_fftgr(:,:,6) 700 if(moved_atm_inside==1)then 701 f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new+ gamma*f_atm(:,:,6) 702 xred(:,:)=f_atm(:,:,1)+ lambda_predict*f_atm(:,:,6) 703 end if 704 705 ! If end_linmin==2, then this vtrial is the good one 706 707 if(end_linmin==2)then 708 709 lambda_old=0.0_dp 710 lambda_new=lambda_predict 711 712 else if(end_linmin==4)then 713 714 ! predict the result of the computation at the trial potential 715 ! defined in the end_linmin==2 case 716 gamma=lambda_predict2/d_lambda_old2 717 ratio=lambda_predict/lambda_new 718 719 ! Take care of vtrial 720 f_fftgr(:,:,1)=vtrial(:,:) 721 722 ABI_ALLOCATE(tmp_fft1,(cplex*nfft,nspden)) 723 ! Take care of vresid 724 tmp_fft1(:,:)=f_fftgr(:,:,i_vresid(2)) 725 f_fftgr(:,:,i_vresid(2))=tmp_fft1(:,:)& 726 & +ratio*(f_fftgr(:,:,i_vresid(1))-tmp_fft1(:,:))& 727 & +gamma*(f_fftgr(:,:,i_vresid(3))-tmp_fft1(:,:)) 728 f_fftgr(:,:,i_vresid(3))=tmp_fft1(:,:) 729 730 ! Take care of rhor 731 tmp_fft1(:,:)=f_fftgr(:,:,i_rhor(2)) 732 f_fftgr(:,:,i_rhor(2))=tmp_fft1(:,:)& 733 & +ratio*(rhor(:,:)-tmp_fft1(:,:))& 734 & +gamma*(f_fftgr(:,:,i_rhor(3))-tmp_fft1(:,:)) 735 f_fftgr(:,:,i_rhor(3))=tmp_fft1(:,:) 736 737 ! Take care of vrespc 738 tmp_fft1(:,:)=f_fftgr(:,:,i_vrespc(2)) 739 f_fftgr(:,:,i_vrespc(2))=tmp_fft1(:,:)& 740 & +ratio*(f_fftgr(:,:,i_vrespc(1))-tmp_fft1(:,:))& 741 & +gamma*(f_fftgr(:,:,i_vrespc(3))-tmp_fft1(:,:)) 742 f_fftgr(:,:,i_vrespc(3))=tmp_fft1(:,:) 743 ABI_DEALLOCATE(tmp_fft1) 744 745 if(moved_atm_inside==1)then 746 do idir=1,3 747 do iatom=1,natom 748 749 ! Take care of xred 750 f_atm(idir,iatom,1)=xred(idir,iatom) 751 752 ! Take care of -HF forces 753 temp=f_atm(idir,iatom,i_vresid(2)) 754 f_atm(idir,iatom,i_vresid(2))=f_atm(idir,iatom,i_vresid(2))& 755 & +ratio*(f_atm(idir,iatom,i_vresid(1))-f_atm(idir,iatom,i_vresid(2)))& 756 & +gamma*(f_atm(idir,iatom,i_vresid(3))-f_atm(idir,iatom,i_vresid(2))) 757 f_atm(idir,iatom,i_vresid(3))=temp 758 759 ! Take care of old xreds 760 temp=f_atm(idir,iatom,i_rhor(2)) 761 f_atm(idir,iatom,i_rhor(2))=f_atm(idir,iatom,i_rhor(2))& 762 & +ratio*( xred(idir,iatom) -f_atm(idir,iatom,i_rhor(2)))& 763 & +gamma*(f_atm(idir,iatom,i_rhor(3))-f_atm(idir,iatom,i_rhor(2))) 764 f_atm(idir,iatom,i_rhor(3))=temp 765 766 ! Take care of preconditioned changes of atomic positions 767 temp=f_atm(idir,iatom,i_vrespc(2)) 768 f_atm(idir,iatom,i_vrespc(2))=f_atm(idir,iatom,i_vrespc(2))& 769 & +ratio*(f_atm(idir,iatom,i_vrespc(1))-f_atm(idir,iatom,i_vrespc(2)))& 770 & +gamma*(f_atm(idir,iatom,i_vrespc(3))-f_atm(idir,iatom,i_vrespc(2))) 771 f_atm(idir,iatom,i_vrespc(3))=temp 772 773 end do 774 end do 775 end if 776 777 ! Since we are at the 2D minimum, the derivative is supposed 778 ! to vanish. Note that dedv_old should not change, by contrast. 779 dedv_old2=0.0_dp 780 d_lambda_old2=-lambda_predict 781 d2edv2_old2=-dedv_old/lambda_predict 782 lambda_old=lambda_predict 783 ilinmin=0 784 785 ! So, jump to the next line 786 iline_cge=iline_cge+1 787 write(message,*)' energy CG update : after 2D interpolation,' 788 call wrtout(std_out,message,'COLL') 789 write(message,*)' computation in the next plane ' 790 call wrtout(std_out,message,'COLL') 791 write(message,*) 792 call wrtout(std_out,message,'COLL') 793 lambda_old=0.0_dp 794 lambda_new=lambda_adapt 795 796 vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2)) 797 if(moved_atm_inside==1)then 798 xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2)) 799 end if 800 801 ! The new trial potential is now generated 802 803 ! End the specific treatment of region IV 804 end if 805 ! 806 ! End the choice between treatment of region I, II, or IV 807 end if 808 809 ! End of choice between initialisation or more developed parts of the CG algorithm 810 else 811 errid = AB7_ERROR_MIXING_ARG 812 errmess = 'scfcge : BUG You should not be here ! ' 813 return 814 end if 815 816 !-------------------------------------- 817 818 !Write information : it will be easy to read by typing grep scfcge logfile 819 820 if(istep==1)then 821 write(message,'(a,a,a)') ' scfcge:',ch10,' scfcge:istep-iline_cge-ilinmin lambda etot resid ' 822 call wrtout(std_out,message,'COLL') 823 end if 824 825 if(ilinmin_input/=0 .or. istep==1)then 826 ! Usual line minimisation step 827 828 if(iline_cge_input<10)then 829 write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' )& 830 & ' scfcge: actual ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input 831 else 832 write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' )& 833 & ' scfcge: actual ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input 834 end if 835 call wrtout(std_out,message,'COLL') 836 837 if( (end_linmin==1.or.end_linmin==-1) .and. istep/=1 )then 838 839 if(end_linmin==1)then 840 write(message, '(a,es13.4,a,i2,a,a)' )& 841 & ' scfcge: predict ',lambda_predict,& 842 & ' suff. close => next line, ilinear=',ilinear,ch10,& 843 & ' scfcge:' 844 else if(end_linmin==-1)then 845 write(message, '(a,es13.4,a,a,a)' )& 846 & ' scfcge: predict ',lambda_predict,& 847 & ' restart the algorithm ',ch10,& 848 & ' scfcge:' 849 end if 850 call wrtout(std_out,message,'COLL') 851 852 if(iline_cge_input<9)then 853 write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' ) & 854 & ' scfcge: start ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old 855 else 856 write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' ) & 857 & ' scfcge: start ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old 858 end if 859 call wrtout(std_out,message,'COLL') 860 861 else if(istep/=1) then 862 write(message, '(a,es13.4,a)' )& 863 & ' scfcge: predict ',lambda_predict,& 864 & ' not close enough => continue minim.' 865 call wrtout(std_out,message,'COLL') 866 end if 867 868 else 869 ! CG prediction 870 if(iline_cge_input<10)then 871 write(message, '(a,i4,a,i1,a,es11.4,es20.12,es12.4,a,i1)' )& 872 & ' scfcge: actual ',istep,'-',iline_cge_input,'-off',& 873 & lambda_adapt,etotal_input,resid_input,', end=',end_linmin 874 else 875 write(message, '(a,i3,a,i2,a,es11.4,es20.12,es12.4,a,i1)' )& 876 & ' scfcge: actual ',istep,'-',iline_cge_input,'-off',& 877 & lambda_adapt,etotal_input,resid_input,', end=',end_linmin 878 end if 879 call wrtout(std_out,message,'COLL') 880 881 if(end_linmin==4)then 882 write(message, '(a)' ) ' scfcge:' 883 call wrtout(std_out,message,'COLL') 884 end if 885 886 end if 887 888 end subroutine scfcge