Statistiques
| Révision :

root / src / egrad_LEPS.f90 @ 8

Historique | Voir | Annoter | Télécharger (3,82 ko)

1
  SUBROUTINE egrad_LEPS(nat,e,geom,grad) 
2
! This program computes the ernergy and gradien in cartesian coordinates
3
! for the cartesian geometry Geom
4

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

    
9
! Number of atoms
10
       INTEGER(KINT), INTENT(IN) :: Nat
11
! Input geometry
12
       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)
13
! output energy and gradient
14
       REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)
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.001d0
21
      
22
! Variables
23
       INTEGER(KINT) :: i,iat,jat
24
       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)
25
 
26
       LOGICAL :: Debug
27
  
28
      REAL(KREAL), external :: ELEPS_xyz
29

    
30

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

    
39
  debug=valid('egrad_leps')
40

    
41
 if (debug) WRITE(*,*) '================ Entering Egrad_leps ==================='
42
 if (debug) THEN
43
    WRITE(*,*) "Cartesian Geometry"
44
    DO I=1,Nat
45
       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3)
46
    END DO
47
 END IF
48
      
49
      ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3))
50

    
51
      e=ELEPS_xyz(nat,Geom)
52

    
53
! We now calculate the gradients using numerical derivatives
54
      Grad=0.d0
55
      GradTmp=0.d0
56
      do iat=1,3
57
         do jat=1,nat
58
            xyztmp=geom
59
            do i=1,6
60
               xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh
61
               gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*ELEPS_xyz(nat,xyztmp)
62
            end do
63
         end do
64
      end do
65
      gradTmp=gradTmp/(60.*hh)
66
      do iat=1,nat
67
       grad(3*iat-2:3*iat)=gradTmp(iat,1:3)
68
      end do
69

    
70
 if (debug) THEN
71
    WRITE(*,*) "Cartesian gradient "
72
    DO I=1,Nat
73
       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3)
74
    END DO
75
 END IF
76

    
77

    
78
      deallocate(xyztmp,gradTmp)
79

    
80

    
81
 if (debug) WRITE(*,*) '================ Exiting Egrad_leps ==================='
82

    
83
! ======================================================================
84
      end 
85

    
86

    
87

    
88
      function ELEPS_xyz(natoms,Xyz)
89

    
90
        use Path_module, only : order
91
        use Io_module, only : au2ev
92

    
93
   IMPLICIT NONE
94
     integer, parameter :: KINT = kind(1)
95
     integer, parameter :: KREAL = kind(1.0d0)
96

    
97

    
98
      INTEGER(KINT) ,INTENT(IN) :: natoms
99
      REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
100
      REAL(KREAL) :: ELEPS_xyz
101

    
102
      INTEGER(KINT) :: i
103
      REAL(KREAL) :: rAB, rBC, rAC
104

    
105
      REAL(KREAL), PARAMETER :: a=0.05, b=0.30, c=0.05
106
      REAL(KREAL), PARAMETER :: dAB=4.746, dBC=4.747, dAC=3.445
107
      REAL(KREAL), PARAMETER :: r0=0.742, alpha=1.942
108

    
109
      REAL(KREAL) :: Qbc,QAc,Qab,Jab,Jbc,Jac
110

    
111
      rAB=0.
112
      rAC=0.
113
      rBC=0.
114
      do i=1,3
115
         rBC=rBC+(xyz(Order(3),i)-xyz(order(2),i))**2
116
         rAB=rAB+(xyz(order(1),i)-xyz(order(2),i))**2
117
         rAC=rAC+(xyz(order(1),i)-xyz(order(3),i))**2
118
      end do
119
      rAB=sqrt(rAB)
120
      rBC=sqrt(rBC)
121
      rAC=sqrt(rAC)
122

    
123
! For consistency with previous articles : rAC is constrained to be 3.742
124
      rAC=3.742
125

    
126
      Qab = dab*(1.5*exp(-2*alpha*(rAB-r0))-exp(-alpha*(rAB-r0)))/(2+2*a)
127
      Qbc = dbc*(1.5*exp(-2*alpha*(rBC-r0))-exp(-alpha*(rBC-r0)))/(2+2*b)
128
      Qac = dac*(1.5*exp (-2*alpha*(rAC-r0))-exp(-alpha*(rAC-r0)))/(2+2*c)
129
      Jab = dab*(exp(-2*alpha*(rAB-r0))-6*exp(-alpha*(rAB-r0)))/(4+4*a)
130
      Jbc = dbc*(exp(-2*alpha*(rBC-r0))-6*exp(-alpha*(rBC-r0)))/(4+4*b)
131
      Jac = dac*(exp(-2*alpha*(rAC-r0))-6*exp(-alpha*(rAC-r0)))/(4+4*c)
132

    
133
! V(x,y) = Qab(x)+Qbc(y)+Qac(rac)-sqrt(Jab(x)**2+Jbc(y)**2+Jac(x)**2-Jab(x)*Jbc(y)-Jbc(y)*Jac(rac)-Jab(x)*Jac(rac))
134

    
135

    
136
      ELEPS_xyz=(Qab+Qbc+Qac-sqrt(Jab**2+Jbc**2+Jac**2-Jab*Jbc-Jbc*Jac-Jab*Jac))/au2eV
137

    
138
      return
139
    end function ELEPS_xyz