TABLE OF CONTENTS
ABINIT/prtxf [ Functions ]
NAME
prtxf
FUNCTION
Compute and print out dimensional cartesian coordinates and forces. Note: for x=cartesian coordinates, t=reduced coordinates (xred), => $ x= R t $ => $ x(1)=rprimd(1,1) t(1)+rprimd(2,1) t(2)+rprimd(3,1) t(3)$ etc. Also $ t = (R^{-1}) x$ . To convert gradients, $d(E)/dx(n) = [d(E)/dt(m)] [dt(m)/dx(n)]$ and $ dt(m)/dx(n) = (R^{-1})_{mn} = G_{nm}$ because G is the inverse transpose of R. Finally then $d(E)/dx(n) = G_{nm} [d(E)/dt(m)]$. The vector $d(E)/dt(m)$ for each atom is input in fred (grad. wrt xred).
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, 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
fred(3,natom)=gradients of Etot (hartree) wrt xred(3,natom) iatfix(3,natom)=1 for each fixed atom along specified direction, else 0 iout=unit number for output file iwfrc=controls force output: 0=> no forces output, 1=>forces out in eV/A and Ha/bohr, 2=>forces out in Ha/bohr natom=number of atoms in unit cell rprimd(3,3)=dimensional real space primitive translations (bohr) xred(3,natom)=relative coordinates of atoms (in terms of prim. transl.)
OUTPUT
(data written to unit iout)
PARENTS
clnup1
CHILDREN
matr3inv,wrtout
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 60 subroutine prtxf(fred,iatfix,iout,iwfrc,natom,rprimd,xred) 61 62 use defs_basis 63 use m_profiling_abi 64 use m_errors 65 66 !This section has been created automatically by the script Abilint (TD). 67 !Do not modify the following lines by hand. 68 #undef ABI_FUNC 69 #define ABI_FUNC 'prtxf' 70 use interfaces_14_hidewrite 71 use interfaces_32_util 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 !scalars 78 integer,intent(in) :: iout,iwfrc,natom 79 !arrays 80 integer,intent(in) :: iatfix(3,natom) 81 real(dp),intent(in) :: fred(3,natom),rprimd(3,3),xred(3,natom) 82 83 !Local variables------------------------------- 84 !scalars 85 integer :: iatom,mu,unfixd 86 real(dp) :: convt,fmax,frms 87 character(len=15) :: format_line21 88 character(len=15) :: format_line25 89 character(len=15) :: format_line 90 character(len=500) :: message 91 !arrays 92 real(dp) :: favg(3),favg_out(3),ff(3),gprimd(3,3),xx(3) 93 94 ! **************************************************************** 95 96 format_line21='(i5,1x,3f21.14)' 97 format_line25='(i5,1x,3f25.14)' 98 99 !Write cartesian coordinates in angstroms 100 call wrtout(iout,' cartesian coordinates (angstrom) at end:','COLL') 101 do iatom=1,natom 102 format_line=format_line21 103 do mu=1,3 104 xx(mu)=(rprimd(mu,1)*xred(1,iatom)+& 105 & rprimd(mu,2)*xred(2,iatom)+& 106 & rprimd(mu,3)*xred(3,iatom))*Bohr_Ang 107 if(xx(mu)>99999 .or. xx(mu)<-9999)format_line=format_line25 108 end do 109 write(message,format_line) iatom,xx 110 call wrtout(iout,message,'COLL') 111 end do 112 113 !Optionally write cartesian forces in eV/Angstrom (also provide same in hartree/bohr) 114 if (iwfrc/=0) then 115 ! First, provide results in hartree/bohr 116 write(message, '(a,a)' ) ch10,' cartesian forces (hartree/bohr) at end:' 117 call wrtout(iout,message,'COLL') 118 frms=zero 119 fmax=zero 120 favg(1)=zero 121 favg(2)=zero 122 favg(3)=zero 123 ! To get cartesian forces from input gradients with respect to 124 ! dimensionless coordinates xred, multiply by G and negate 125 ! (see notes at top of this subroutine) 126 call matr3inv(rprimd,gprimd) 127 ! First compute (spurious) average force favg 128 do iatom=1,natom 129 do mu=1,3 130 ff(mu)=-(gprimd(mu,1)*fred(1,iatom)+& 131 & gprimd(mu,2)*fred(2,iatom)+& 132 & gprimd(mu,3)*fred(3,iatom)) 133 favg(mu)=favg(mu)+ff(mu) 134 end do 135 end do 136 favg(1) = favg(1)/dble(natom) 137 favg(2) = favg(2)/dble(natom) 138 favg(3) = favg(3)/dble(natom) 139 140 ! Subtract off average force in what follows 141 ! (avg is also subtracted off in carfor, called by loopcv, 142 ! called by grad) 143 unfixd=0 144 do iatom=1,natom 145 format_line=format_line21 146 do mu=1,3 147 ff(mu)=-(gprimd(mu,1)*fred(1,iatom)+& 148 & gprimd(mu,2)*fred(2,iatom)+& 149 & gprimd(mu,3)*fred(3,iatom))-favg(mu) 150 if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25 151 ! For rms and max force, include only unfixed components 152 if (iatfix(mu,iatom) /= 1) then 153 unfixd=unfixd+1 154 frms=frms+ff(mu)**2 155 fmax=max(fmax,abs(ff(mu))) 156 end if 157 end do 158 write(message, format_line) iatom,ff 159 call wrtout(iout,message,'COLL') 160 end do 161 if ( unfixd /= 0 ) frms = sqrt(frms/dble(unfixd)) 162 163 ! The average force is obtained from the cancellation of numbers 164 ! of typical size unity, so an absolute value lower 165 ! than tol14 is meaningless for the output file. 166 favg_out(:)=favg(:) 167 if(abs(favg_out(1))<tol14)favg_out(1)=zero 168 if(abs(favg_out(2))<tol14)favg_out(2)=zero 169 if(abs(favg_out(3))<tol14)favg_out(3)=zero 170 171 write(message, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',frms,fmax,favg_out(1:3),' h/b' 172 call wrtout(iout,message,'COLL') 173 174 if (iwfrc==1) then 175 176 write(message, '(a,a)' )ch10,' cartesian forces (eV/Angstrom) at end:' 177 call wrtout(iout,message,'COLL') 178 convt=Ha_eV/Bohr_Ang 179 180 ! Note: subtract off average force 181 do iatom=1,natom 182 format_line=format_line21 183 do mu=1,3 184 ff(mu)=(-(gprimd(mu,1)*fred(1,iatom)+& 185 & gprimd(mu,2)*fred(2,iatom)+& 186 & gprimd(mu,3)*fred(3,iatom))-favg(mu))*convt 187 if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25 188 end do 189 write(message, format_line) iatom,ff 190 call wrtout(iout,message,'COLL') 191 end do 192 write(message, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',convt*frms,convt*fmax,convt*favg_out(1:3),' e/A' 193 call wrtout(iout,message,'COLL') 194 195 end if 196 end if 197 198 end subroutine prtxf