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