1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8subroutine ofc(fa, dx, na, has_constr, first) 9! ******************************************************************* 10! Writes force constants matrix to file 11! Input forces are in Ry/Bohr and input displacements are in Bohr. 12! Force Constants written in file are in eV / Ang ** 2 13! Written by P.Ordejon. August'98. 14! Dynamic memory and save attribute for fres introduced by J.Gale 15! Sept'99. 16! Re-structured by Alberto Garcia, April 2007 17! Added constrained by Nick Papior, June 2015 18! Added displacement value by Nick Papior, June 2018 19! ********* INPUT *************************************************** 20! real*8 fa(3,na) : atomic forces (in Ry / Bohr) 21! real*8 dx : atomic displacements (in Bohr) 22! integer na : number of atoms 23! logical has_constr : whether the forces are constrained or not 24! logical first : first call (initializes the force) 25! ********** BEHAVIOUR ********************************************** 26! On the first call (undisplaced coordinates), the forces should be 27! zero (relaxed structure). 28! However, since the relaxation is usually not perfect, some 29! residual forces are obtained. These residual forces are 30! substracted from the forces on other steps, to calculate the 31! force constants matrix 32! ******************************************************************* 33 34 use precision 35 use files, only : slabel, label_length 36 use units, only : Ang, eV 37 use alloc, only : re_alloc 38 39 implicit none 40 41 integer, intent(in) :: na 42 real(dp), intent(in) :: dx, fa(3,na) 43 logical, intent(in) :: has_constr, first 44 45 external :: io_assign, io_close 46 47 ! Local reference force for the un-displaced configuration. 48 real(dp), dimension(:,:,:), pointer, save :: fres => null() 49 50 ! Saved variables and arrays 51 character(len=label_length+4) :: fname 52 real(dp), dimension(:,:), pointer :: fr 53 real(dp) :: tmp 54 55 integer :: i, ix, iu 56 57 if ( .not. associated(fres) ) then 58 call re_alloc( fres, 1, 3, 1, na, 1, 2, name='fres', routine='ofc' ) 59 end if 60 61 if ( has_constr ) then 62 fname = trim(slabel) // '.FCC' 63 fr => fres(:,:,2) 64 else 65 fname = trim(slabel) // '.FC' 66 fr => fres(:,:,1) 67 end if 68 69 if ( first ) then 70 71 call io_assign(iu) 72 open( iu, file=fname, status='unknown' ) 73 rewind(iu) 74 75 tmp = dx / Ang 76 if ( has_constr ) then 77 write(iu,'(2a,i0,tr1,e22.16)') 'Force constants matrix (constrained). ', & 78 'n_atoms, displacement [Ang]: ', na, tmp 79 else 80 write(iu,'(2a,i0,tr1,e22.16)') 'Force constants matrix. ', & 81 'n_atoms, displacement [Ang]: ', na, tmp 82 end if 83 84 call io_close(iu) 85 86 ! Copy over the residual (zero point) forces 87 fr(:,:) = fa(:,:) 88 89 ! We should not save anything now 90 return 91 92 end if 93 94 call io_assign(iu) 95 open( iu, file=fname, status='old',position="append", action="write") 96 97 tmp = Ang ** 2 / eV / dx 98 do i = 1 , na 99 write(iu,'(3(tr1,e17.9))') ((-fa(ix,i)+fr(ix,i))*tmp,ix=1,3) 100 enddo 101 102 call io_close(iu) 103 104end subroutine ofc 105