Statistiques
| Révision :

root / src / egrad_test_2D.f90 @ 8

Historique | Voir | Annoter | Télécharger (4,13 ko)

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 Io_module, only : IOTMP
130

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

    
135

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

    
141
      INTEGER(KINT) :: i,j
142

    
143
      LOGICAL :: Debug
144
      LOGICAL, SAVE :: First=.TRUE.
145
 
146
  INTERFACE
147
     function valid(string) result (isValid)
148
       CHARACTER(*), intent(in) :: string
149
       logical                  :: isValid
150
     END function VALID
151
     
152
  END INTERFACE
153

    
154
  debug=valid('e2D')
155

    
156
 if (First) THEN
157
! we read the data
158
   open(IOTMP,File="CCSDT.subtr.reg")
159
   DO I=1, 501
160
    DO j=1,501
161
     READ(IOTMP,*) x,y,E
162
     Data(int(x*100+0.1)+1,int(y*100+0.1)+1)=E
163
     write(*,*) x,y,int(x*100+0.1)+1, int(y*100+0.1)+1,E
164
    END DO
165
! we skip a blank line
166
     READ(IOTMP,*) 
167
  END DO
168
  CLOSE(IOTMP)
169
  First=.FALSE.
170
 END IF
171

    
172
    x=Xyz(2,1)
173
    y=Xyz(3,1)
174
    E2D=Data(int(x*100+0.1)+1,int(y*100+0.1)+1)
175

    
176
      return
177
    end function E2D