TABLE OF CONTENTS
ABINIT/constrf [ Functions ]
NAME
constrf
FUNCTION
Computes projected forces, fpcart, which satisfy a set of constraint equations of the form Sum[mu,iatom]: wtatcon(mu,iatom,iconeq)*fpcart(mu,iatom) = 0 (iconeq=1,nconeq). These projected forces are returned in fcart and thus replace the original forces.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (SCE, 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
iatfix(3,natom)=1 for frozen atom along each direction, 0 for unfrozen natom=number of atoms in cell nconeq=number of atomic constraint equations prtvol=control print volume and debugging rprimd(3,3)=dimensional primitive translations in real space (bohr) wtatcon(3,natom,nconeq)=weights for atomic constraints xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
diffor=maximum absolute change in component of projected fcart between present
and previous SCF cycle
fred(3,natom)=grads of Etot wrt reduced coordinates (hartree)
maxfor=maximum absolute value of fcart
SIDE EFFECTS
fcart(3,natom)=cartesian forces (hartree/bohr) on input, projected forces on output forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr)
TODO
PARENTS
forces
CHILDREN
dposv,prtxvf,wrtout,xred2xcart
SOURCE
49 #if defined HAVE_CONFIG_H 50 #include "config.h" 51 #endif 52 53 #include "abi_common.h" 54 55 subroutine constrf(diffor,fcart,forold,fred,iatfix,ionmov,maxfor,natom,& 56 & nconeq,prtvol,rprimd,wtatcon,xred) 57 58 use defs_basis 59 use m_linalg_interfaces 60 use m_errors 61 use m_profiling_abi 62 63 use m_geometry, only : xred2xcart 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'constrf' 69 use interfaces_14_hidewrite 70 use interfaces_67_common, except_this_one => constrf 71 !End of the abilint section 72 73 implicit none 74 75 !Arguments ------------------------------------ 76 !scalars 77 integer,intent(in) :: ionmov,natom,nconeq,prtvol 78 real(dp),intent(out) :: diffor,maxfor 79 !arrays 80 integer,intent(in) :: iatfix(3,natom) 81 real(dp),intent(in) :: rprimd(3,3),wtatcon(3,natom,nconeq) 82 real(dp),intent(inout) :: fcart(3,natom),forold(3,natom),xred(3,natom) 83 real(dp),intent(inout) :: fred(3,natom) !vz_i 84 85 !Local variables ------------------------- 86 !scalars 87 integer :: iatom,iconeq,iconeq1,iconeq2,index,info,mu,prtvel 88 character(len=500) :: message 89 !arrays 90 real(dp),allocatable :: fcartvec(:),fpcart(:,:),fvector(:),vel_dummy(:,:) 91 real(dp),allocatable :: wmatrix(:,:),wtatconvec(:,:),wtcoeffs(:),xcart(:,:) 92 93 !************************************************************************ 94 95 !Allocate temporary variables 96 ABI_ALLOCATE(fpcart,(3,natom)) 97 ABI_ALLOCATE(fcartvec,(3*natom)) 98 ABI_ALLOCATE(fvector,(nconeq)) 99 ABI_ALLOCATE(vel_dummy,(3,natom)) 100 ABI_ALLOCATE(wmatrix,(nconeq,nconeq)) 101 ABI_ALLOCATE(wtatconvec,(3*natom,nconeq)) 102 ABI_ALLOCATE(wtcoeffs,(nconeq)) 103 ABI_ALLOCATE(xcart,(3,natom)) 104 105 !If prtvol>10, output coordinates and forces prior to projecting 106 if(prtvol>=10)then 107 write(message,'(a)')' constrf - coordinates and forces prior to constraint projections:' 108 call wrtout(std_out,message,'COLL') 109 call xred2xcart(natom,rprimd,xcart,xred) 110 prtvel=0 111 call prtxvf(fcart,fred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred) 112 end if 113 114 !Transfer fcart and wtatcon to flat vectors 115 index=0 116 do iatom=1,natom 117 do mu=1,3 118 index=index+1 119 fcartvec(index)=fcart(mu,iatom) 120 wtatconvec(index,:)=wtatcon(mu,iatom,:) 121 end do 122 end do 123 124 !Compute a matrix (wmatrix) and vector (fvector) such that solving 125 !the linear equations wmatrix*wcoeffs=fvector gives the coefficients 126 !of wtatcon (wcoeffs) needed to compute the projected forces 127 do iconeq2=1,nconeq 128 fvector(iconeq2)=ddot(3*natom,fcartvec,1,wtatconvec(1,iconeq2),1) 129 do iconeq1=1,nconeq 130 wmatrix(iconeq1,iconeq2)=ddot(3*natom,wtatconvec(1,iconeq1),1,wtatconvec(1,iconeq2),1) 131 end do 132 end do 133 134 !Solve the system of linear equations, wmatrix*wcoeffs=fvector 135 call dposv('U',nconeq,1,wmatrix,nconeq,fvector,nconeq,info) 136 137 if (info/=0) then 138 write(message, '(a,a,a,a,a)' )& 139 & ' Constraint matrix is not positive definite,',ch10,& 140 & ' probably because constraints are linearly dependent.',ch10,& 141 & ' Action : Check for linear dependence of constraints.' 142 MSG_ERROR(message) 143 end if 144 145 !The solution vector is returned in fvector, so copy it to a more sensible location 146 wtcoeffs(:)=fvector(:) 147 148 !Compute the projected forces, which now satisfy all the constraint equations 149 fpcart(:,:)=fcart(:,:) 150 do iconeq=1,nconeq 151 fpcart(:,:)=fpcart(:,:)-wtcoeffs(iconeq)*wtatcon(:,:,iconeq) 152 end do 153 154 !Reconvert constrained forces back from fpcart to fred 155 do iatom=1,natom 156 do mu=1,3 157 fred(mu,iatom)= - (rprimd(1,mu)*fpcart(1,iatom)+& 158 & rprimd(2,mu)*fpcart(2,iatom)+& 159 & rprimd(3,mu)*fpcart(3,iatom)) 160 end do 161 end do 162 163 !If prtvol>=10, output coordinates and forces after projecting 164 if(prtvol>=10)then 165 write(message,'(a)')' constrf - coordinates and forces after constraint projections:' 166 call wrtout(std_out,message,'COLL') 167 prtvel=0 168 call prtxvf(fpcart,fred,iatfix,06,natom,prtvel,vel_dummy,xcart,xred) 169 end if 170 171 !Copy the constrained forces, fpcart, back to fcart 172 fcart(:,:)=fpcart(:,:) 173 174 !Compute maximal force and maximal difference of the projected forces, 175 !overriding the values already computed in forces 176 maxfor=0.0_dp 177 diffor=0.0_dp 178 do iatom=1,natom 179 do mu=1,3 180 if (iatfix(mu,iatom) /= 1) then 181 maxfor=max(maxfor,abs(fcart(mu,iatom))) 182 diffor=max(diffor,abs(fcart(mu,iatom)-forold(mu,iatom))) 183 else if (ionmov==4 .or. ionmov==5) then 184 ! Make the force vanish on fixed atoms when ionmov=4 or 5 185 ! This is because fixing of atom cannot be imposed at the 186 ! level of a routine similar to brdmin or moldyn for these options. 187 fcart(mu,iatom)=0.0_dp 188 end if 189 end do 190 end do 191 forold(:,:)=fcart(:,:) 192 193 !Dellocate temporary variables 194 ABI_DEALLOCATE(fpcart) 195 ABI_DEALLOCATE(fcartvec) 196 ABI_DEALLOCATE(fvector) 197 ABI_DEALLOCATE(vel_dummy) 198 ABI_DEALLOCATE(wmatrix) 199 ABI_DEALLOCATE(wtatconvec) 200 ABI_DEALLOCATE(wtcoeffs) 201 ABI_DEALLOCATE(xcart) 202 203 end subroutine constrf