← Back to branch summary

~siesta-sip/siesta/sip-solver

« back to all changes in this revision

Viewing changes to Src/ofc.f

  • Committer: Alberto Garcia
  • Date: 2015-06-22 09:32:02 UTC
  • Revision ID: albertog@icmab.es-20150622093202-enhjztqfku2e976q
Revert previous change affecting FA and FC files. FAC/FCC files.

* At revno 464, a change was introduced to output the constrained
  forces to the FA and FC files.  This had unintended consequences for
  some analysis tools.

  The change has been reverted, and the FC/FA files are now created as
  before. In addition, if constraints are used, new FCC/FAC files
  will be created with the constrained forces.

  The vibra utility will default to use the FCC files if the
  GeometryConstraints block exists, and in any case the
  force-constants file to use can be specified with the Vibra.FC fdf
  option.

  To support these changes, the write_forces subroutine has now 
  the geometry step as an additional argument.

  The writing of forces has been moved to its appropriate place in the
  write_subs modules.

* Some more cosmetic changes in the Vibra package.

modified:
  Docs/siesta.tex
  Src/Makefile
  Src/iofa.f
  Src/ofc.f
  Src/siesta_analysis.F
  Src/siesta_forces.F
  Src/state_analysis.F
  Src/write_subs.F
  Util/Vibra/Docs/CHANGES
  Util/Vibra/Docs/vibra.tex
  Util/Vibra/Examples/README
  Util/Vibra/Examples/si54.bands
  Util/Vibra/Examples/si54.fdf
  Util/Vibra/Src/Makefile
  Util/Vibra/Src/fcbuild.f
  Util/Vibra/Src/recoor.f
  Util/Vibra/Src/vibra.f

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/ofc.f
8
8
! Use of this software constitutes agreement with the full conditions
9
9
! given in the SIESTA license, as signed by all legitimate users.
10
10
!
11
 
      subroutine ofc(fa, dx, na)
 
11
      subroutine ofc(fa, dx, na, has_constr, first)
12
12
C *******************************************************************
13
13
C Writes force constants matrix to file
14
14
C Input forces are in Ry/Bohr and input displacements are in Bohr.
17
17
C Dynamic memory and save attribute for fres introduced by J.Gale
18
18
C Sept'99.
19
19
C Re-structured by Alberto Garcia, April 2007
 
20
C Added constrained by Nick Papior, June 2015    
20
21
C ********* INPUT ***************************************************
21
22
C real*8 fa(3,na)             : atomic forces (in Ry / Bohr)
22
23
C real*8 dx                   : atomic displacements (in Bohr)
23
24
C integer na                  : number of atoms
 
25
c logical has_constr          : whether the forces are constrained or not
 
26
c logical first               : first call (initializes the force)
24
27
C ********** BEHAVIOUR **********************************************
25
28
C On the first call (undisplaced coordinates), the forces should be 
26
29
C zero (relaxed structure).
37
40
 
38
41
      implicit          none
39
42
 
40
 
      integer, intent(in)  ::    na
41
 
      real(dp), intent(in) ::    dx, fa(3,na)
 
43
      integer,  intent(in) :: na
 
44
      real(dp), intent(in) :: dx, fa(3,na)
 
45
      logical,  intent(in) :: has_constr, first
42
46
 
43
47
      external          io_assign, io_close
44
48
 
45
49
C Saved  variables and arrays
46
 
      character(len=label_length+3), save :: fname
47
 
      logical,                       save :: frstme = .true.
48
 
      real(dp), dimension(:,:), pointer, save :: fres
49
 
 
50
 
      integer    :: i, ix, unit1
51
 
 
52
 
 
53
 
      if (frstme) then
54
 
 
55
 
        fname = trim(slabel) // '.FC'
56
 
        nullify( fres )
57
 
        call re_alloc( fres, 1, 3, 1, na, name='fres', routine='ofc' )
58
 
 
59
 
        call io_assign(unit1)
60
 
        open( unit1, file=fname, status='unknown' )
61
 
        rewind(unit1)
62
 
        write(unit1,'(a)') 'Force constants matrix'
63
 
        call io_close(unit1)
64
 
        fres(:,:) = fa(:,:)
65
 
        frstme = .false.
66
 
 
 
50
      character(len=label_length+4) :: fname
 
51
      real(dp), dimension(:,:,:), pointer, save :: fres => null()
 
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,
 
59
     &        name='fres', routine='ofc' )
 
60
      end if
 
61
 
 
62
      if ( has_constr ) then
 
63
         fname = trim(slabel) // '.FCC'
 
64
         fr => fres(:,:,2)
67
65
      else
68
 
 
69
 
         call io_assign(unit1)
70
 
         open( unit1, file=fname, status='old',position="append",
71
 
     $         action="write")
72
 
         do i=1,na
73
 
            write(unit1,'(3f15.7)') ((-fa(ix,i)+fres(ix,i))*
74
 
     .           Ang**2/eV/dx, ix=1,3)
75
 
         enddo
76
 
         call io_close(unit1)
77
 
 
78
 
      endif
79
 
 
 
66
         fname = trim(slabel) // '.FC'
 
67
         fr => fres(:,:,1)
 
68
      end if
 
69
 
 
70
      if ( first ) then
 
71
 
 
72
         call io_assign(iu)
 
73
         open( iu, file=fname, status='unknown' )
 
74
         rewind(iu)
 
75
 
 
76
         if ( has_constr ) then
 
77
            write(iu,'(a)') 'Force constants matrix (constrained)'
 
78
         else
 
79
            write(iu,'(a)') 'Force constants matrix'
 
80
         end if
 
81
 
 
82
         call io_close(iu)
 
83
 
 
84
         ! Copy over the residual (zero point) forces
 
85
         fr(:,:) = fa(:,:)
 
86
 
 
87
         ! We should not save anything now
 
88
         return
 
89
 
 
90
      end if
 
91
 
 
92
      call io_assign(iu)
 
93
      open( iu, file=fname, status='old',position="append",
 
94
     $     action="write")
 
95
 
 
96
      tmp = Ang ** 2 / eV / dx
 
97
      do i = 1 , na
 
98
         write(iu,'(3f15.7)') ((-fa(ix,i)+fr(ix,i))*tmp,ix=1,3)
 
99
      enddo
 
100
      
 
101
      call io_close(iu)
 
102
      
80
103
      end