Statistics
| Revision:

root / src / egrad_test_2D.f90 @ 5

History | View | Annotate | Download (4.4 kB)

1
  SUBROUTINE egrad_test_2D(nat,e,geom,grad) 
2
! This program computes the energy and the gradient
3
! to cheat PATH we use a 3 atoms system
4

    
5
   IMPLICIT NONE
6
     integer, parameter :: KINT = kind(1)
7
     integer, parameter :: KREAL = kind(1.0d0)
8

    
9
       INTEGER(KINT), INTENT(IN) :: Nat
10
       REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)
11
       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)
12

    
13
!     Bohr --> Angstr
14
!      real(KREAL), parameter :: BOHR   = 0.52917726D+00
15
!
16

    
17
! Parameters to define the surface
18
      INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/)
19
      INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/)
20
      REAL(KREAL), PARAMETER :: hh=0.01d0
21
      INTEGER(KINT), PARAMETER :: IOTMP=99
22
! Variables
23
       INTEGER(KINT) :: i, iat, jat,j
24
       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)
25

    
26
       LOGICAL :: Debug
27
       LOGICAL, save :: First
28
  
29
      REAL(KREAL) :: x,y
30

    
31

    
32

    
33
  INTERFACE
34
     function valid(string) result (isValid)
35
       CHARACTER(*), intent(in) :: string
36
       logical                  :: isValid
37
     END function VALID
38

    
39

    
40
      function E2D(natoms,Xyz)
41

    
42
        use Path_module, only : order
43
        use Io_module, only : IOTMP
44

    
45
   IMPLICIT NONE
46
     integer, parameter :: KINT = kind(1)
47
     integer, parameter :: KREAL = kind(1.0d0)
48

    
49

    
50
      INTEGER(KINT) ,INTENT(IN) :: natoms
51
      REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
52
      REAL(KREAL) :: E2D
53
    END function E2D
54

    
55
     
56
  END INTERFACE
57

    
58
  debug=valid('egrad_test_2D')
59

    
60
 if (debug) WRITE(*,*) '================ Entering Egrad_test ==================='
61
 if (debug) THEN
62
    WRITE(*,*) "Cartesian Geometry"
63
    DO I=1,Nat
64
       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3)
65
    END DO
66
 END IF
67
      
68
      ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3))
69

    
70
      if (first) THEN
71
         OPEN(IOTMP,File="tt")
72
         DO I=1, 501
73
            DO j=1,501
74
               x=(i-1.)/100.
75
               y=(j-1.)/100.
76
               xyztmp(2,1)=x
77
               xyztmp(3,1)=y
78
               write(IOTMP,*) x,y,E2D(nat,xyztmp)
79
            END DO
80
! we skip a blank line
81
            WRITe(IOTMP,*) 
82
         END DO
83
         CLOSE(IOTMP)
84
         First=.FALSE.
85
      END IF
86

    
87
      e=E2D(nat,Geom)
88

    
89
! We now calculate the gradients using numerical derivatives
90
      Grad=0.d0
91
      GradTmp=0.d0
92
      do iat=1,1
93
         do jat=2,nat
94
            xyztmp=geom
95
            do i=1,6
96
               xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh
97
               gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*E2D(nat,xyztmp)
98
            end do
99
         end do
100
      end do
101
      gradTmp=gradTmp/(60.*hh)
102
      do iat=1,nat
103
       grad(3*iat-2:3*iat)=gradTmp(iat,1:3)
104
      end do
105

    
106
 if (debug) THEN
107
    WRITE(*,*) "Cartesian gradient GradTmp"
108
    DO I=1,Nat
109
       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3)
110
    END DO
111
    WRITE(*,*) "Cartesian gradient Grad"
112
    WRITE(*,'(50(1X,F12.6))') Grad
113

    
114
 END IF
115

    
116

    
117
      deallocate(xyztmp,gradTmp)
118

    
119

    
120
 if (debug) WRITE(*,*) '================ Exiting Egrad_test_2D ==================='
121

    
122
! ======================================================================
123
   end SUBROUTINE egrad_test_2D
124

    
125

    
126

    
127
      function E2D(natoms,Xyz)
128

    
129
        use Path_module, only : order
130
        use Io_module, only : IOTMP
131

    
132
   IMPLICIT NONE
133
     integer, parameter :: KINT = kind(1)
134
     integer, parameter :: KREAL = kind(1.0d0)
135

    
136

    
137
      INTEGER(KINT) ,INTENT(IN) :: natoms
138
      REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
139
      REAL(KREAL) :: E,x,y, E2D
140
      REAL(KREAL), save :: Data(501,501)
141

    
142
      INTEGER(KINT) :: i,j
143
      REAL(KREAL) :: dCH,dNH,dCN
144

    
145
!      real(KREAL) :: autoA,autoeV
146
      real(KREAL) :: autoeV
147
!      parameter (autoA=0.52917715d0)
148
      parameter (autoeV=27.21183d0)
149

    
150

    
151
      real(KREAL) :: r1, r2
152

    
153
      LOGICAL :: Debug
154
      LOGICAL, SAVE :: First=.TRUE.
155
 
156
  INTERFACE
157
     function valid(string) result (isValid)
158
       CHARACTER(*), intent(in) :: string
159
       logical                  :: isValid
160
     END function VALID
161
     
162
  END INTERFACE
163

    
164
  debug=valid('e2D')
165

    
166
 if (First) THEN
167
! we read the data
168
   open(IOTMP,File="CCSDT.subtr.reg")
169
   DO I=1, 501
170
    DO j=1,501
171
     READ(IOTMP,*) x,y,E
172
     Data(int(x*100+0.1)+1,int(y*100+0.1)+1)=E
173
     write(*,*) x,y,int(x*100+0.1)+1, int(y*100+0.1)+1,E
174
    END DO
175
! we skip a blank line
176
     READ(IOTMP,*) 
177
  END DO
178
  CLOSE(IOTMP)
179
  First=.FALSE.
180
 END IF
181

    
182
    x=Xyz(2,1)
183
    y=Xyz(3,1)
184
    E2D=Data(int(x*100+0.1)+1,int(y*100+0.1)+1)
185

    
186
      return
187
    end function E2D