TABLE OF CONTENTS
ABINIT/findmin [ Functions ]
NAME
findmin
FUNCTION
Compute the minimum of a function whose value
and derivative are known at two points.
Also deduce different quantities at this predicted
point, and at the two other points
It uses a quartic interpolation, with the supplementary
condition that the second derivative vanishes at one and
only one point (See Schlegel, J. Comp. Chem. 3, 214 (1982).
For this option, lambda_1 must be 1 (new point),
and lambda_2 must be 0 (old point).
Also, if the derivative at the new point is more negative
than the derivative at the old point, the predicted
point cannot correspond to a minimum, but will be lambda=2.5_dp,
if the energy of the second point is lower than the energy
of the first point.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (XG, GMR) 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
etotal_1=first value of the function etotal_2=second value of the function dedv_1=first value of the derivative dedv_2=second value of the derivative lambda_1=first value of the argument lambda_2=second value of the argument
OUTPUT
dedv_predict=predicted value of the derivative (usually zero, except if choice=4, if it happens that a minimum cannot be located, and a trial step is taken) d2edv2_predict=predicted value of the second derivative (not if choice=4) d2edv2_1=first value of the second derivative (not if choice=4) d2edv2_2=second value of the second derivative (not if choice=4) etotal_predict=predicted value of the function lambda_predict=predicted value of the argument status= 0 if everything went normally ; 1 if negative second derivative 2 if some other problem
PARENTS
m_bfgs
CHILDREN
wrtout
SOURCE
60 #if defined HAVE_CONFIG_H 61 #include "config.h" 62 #endif 63 64 #include "abi_common.h" 65 66 67 subroutine findmin(dedv_1,dedv_2,dedv_predict,& 68 & d2edv2_1,d2edv2_2,d2edv2_predict,& 69 & etotal_1,etotal_2,etotal_predict,& 70 & lambda_1,lambda_2,lambda_predict,status) 71 72 use defs_basis 73 use m_errors 74 75 !This section has been created automatically by the script Abilint (TD). 76 !Do not modify the following lines by hand. 77 #undef ABI_FUNC 78 #define ABI_FUNC 'findmin' 79 use interfaces_14_hidewrite 80 !End of the abilint section 81 82 implicit none 83 84 !Arguments ------------------------------------ 85 !scalars 86 integer,intent(out) :: status 87 real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2 88 real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict 89 real(dp),intent(out) :: etotal_predict,lambda_predict 90 91 !Local variables------------------------------- 92 !scalars 93 real(dp) :: aa,bb,bbp,cc,ccp,d_lambda,dd 94 real(dp) :: discr,ee,eep,lambda_shift,sum1,sum2,sum3,uu 95 real(dp) :: uu3,vv,vv3 96 character(len=500) :: message 97 98 ! ************************************************************************* 99 100 !DEBUG 101 !write(std_out,*)' findmin : enter' 102 !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2 103 !ENDDEBUG 104 105 status=0 106 d_lambda=lambda_1-lambda_2 107 108 !DEBUG 109 !do choice=3,1,-1 110 !ENDDEBUG 111 112 if(abs(lambda_1-1.0_dp)>tol12 .or. abs(lambda_2)>tol12) then 113 message = ' For choice=4, lambda_1 must be 1 and lambda_2 must be 0.' 114 MSG_BUG(message) 115 end if 116 117 !Evaluate quartic interpolation 118 !etotal = aa + bb * lambda + cc * lambda**2 + dd * lambda**3 + ee * lambda**4 119 !Impose positive second derivative everywhere, with 120 !one point where it vanishes : 3*dd**2=8*cc*ee 121 aa=etotal_2 122 bb=dedv_2 123 sum1=etotal_1-aa-bb 124 sum2=dedv_1-bb 125 sum3=sum2-2.0_dp*sum1 126 127 !Build the discriminant of the associated 2nd degree equation 128 discr=sum2**2-3.0_dp*sum3**2 129 if(discr<0.0_dp .or. sum2<0.0_dp)then 130 131 ! jmb init 132 d2edv2_2=0.0 133 d2edv2_1=0.0 134 d2edv2_predict=0.0 135 136 ! Even if there is a problem, try to keep going ... 137 message = 'The 2nd degree equation has no positive root (choice=4).' 138 MSG_WARNING(message) 139 status=2 140 if(etotal_1<etotal_2)then 141 write(message, '(a,a,a)' )& 142 & 'Will continue, since the new total energy is lower',ch10,& 143 & 'than the old. Take a larger step in the same direction.' 144 MSG_COMMENT(message) 145 lambda_predict=2.5_dp 146 else 147 write(message, '(a,a,a,a,a)' )& 148 & 'There is a problem, since the new total energy is larger',ch10,& 149 & 'than the old (choice=4).',ch10,& 150 & 'I take a point between the old and new, close to the old .' 151 MSG_COMMENT(message) 152 lambda_predict=0.25_dp 153 end if 154 ! Mimick a zero-gradient lambda, in order to avoid spurious 155 ! action of the inverse hessian (the next line would be a realistic estimation) 156 dedv_predict=0.0_dp 157 ! dedv_predict=dedv_2+lambda_predict*(dedv_1-dedv_2) 158 ! Uses the energies, and the gradient at lambda_2 159 etotal_predict=etotal_2+dedv_2*lambda_predict& 160 & +(etotal_1-etotal_2-dedv_2)*lambda_predict**2 161 162 else 163 164 ! Here, there is an acceptable solution to the 2nd degree equation 165 discr=sqrt(discr) 166 ! The root that gives the smallest ee corresponds to -discr 167 ! This is the one to be used: one aims at modelling the 168 ! behaviour of the function as much as possible with the 169 ! lowest orders of the polynomial, not the quartic term. 170 ee=(sum2-discr)*0.5_dp 171 dd=sum3-2.0_dp*ee 172 cc=sum1-dd-ee 173 174 ! DEBUG 175 ! write(std_out,*)'aa,bb,cc,dd,ee',aa,bb,cc,dd,ee 176 ! ENDDEBUG 177 178 ! Now, must find the unique root of 179 !$0 = bb + 2*cc * lambda + 3*dd * lambda^2 + 4*ee * lambda^3$ 180 ! This root is unique because it was imposed that the second derivative 181 ! of the quartic polynomial is everywhere positive. 182 ! First, remove the quadratic term, by a shift of lambda 183 ! lambdap=lambda-lambda_shift 184 !$0 = bbp + ccp * lambdap + eep * lambdap^3$ 185 eep=4.0_dp*ee 186 lambda_shift=-dd/(4.0_dp*ee) 187 ccp=2.0_dp*cc-12.0_dp*ee*lambda_shift**2 188 bbp=bb+ccp*lambda_shift+eep*lambda_shift**3 189 190 ! DEBUG 191 ! write(std_out,*)'bbp,ccp,eep,lambda_shift',bbp,ccp,eep,lambda_shift 192 ! ENDDEBUG 193 194 ! The solution of a cubic polynomial equation is as follows : 195 discr=(bbp/eep)**2+(4.0_dp/27.0_dp)*(ccp/eep)**3 196 ! In the present case, discr will always be positive 197 discr=sqrt(discr) 198 uu3=0.5_dp*(-bbp/eep+discr) ; uu=sign((abs(uu3))**(1.0_dp/3.0_dp),uu3) 199 vv3=0.5_dp*(-bbp/eep-discr) ; vv=sign((abs(vv3))**(1.0_dp/3.0_dp),vv3) 200 lambda_predict=uu+vv 201 202 ! Restore the shift 203 lambda_predict=lambda_predict+lambda_shift 204 etotal_predict=aa+bb*lambda_predict+cc*lambda_predict**2+& 205 & dd*lambda_predict**3+ee*lambda_predict**4 206 dedv_predict=bb+2.0_dp*cc*lambda_predict+3.0_dp*dd*lambda_predict**2+& 207 & 4.0_dp*ee*lambda_predict**3 208 d2edv2_1=2*cc+6*dd*lambda_1+12*ee*lambda_1**2 209 d2edv2_2=2*cc+6*dd*lambda_2+12*ee*lambda_2**2 210 d2edv2_predict=2*cc+6*dd*lambda_predict+12*ee*lambda_predict**2 211 212 end if 213 214 write(message, '(a,i3)' )' line minimization, algorithm ',4 215 call wrtout(std_out,message,'COLL') 216 write(message, '(a,a)' )' lambda etotal ',' dedv d2edv2 ' 217 call wrtout(std_out,message,'COLL') 218 write(message, '(a,es12.4,es18.10,2es12.4)' )' old point :',lambda_2,etotal_2,dedv_2,d2edv2_2 219 call wrtout(std_out,message,'COLL') 220 write(message, '(a,es12.4,es18.10,2es12.4)' )' new point :',lambda_1,etotal_1,dedv_1,d2edv2_1 221 call wrtout(std_out,message,'COLL') 222 write(message, '(a,es12.4,es18.10,2es12.4)' )' predicted point :',lambda_predict,etotal_predict,dedv_predict,d2edv2_predict 223 call wrtout(std_out,message,'COLL') 224 write(message, '(a)' ) ' ' 225 call wrtout(std_out,message,'COLL') 226 227 end subroutine findmin