Statistiques
| Révision :

root / src / egrad_vasp.f90 @ 8

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

1
subroutine egrad_vasp(e,geomcart,gradcart)
2

    
3
  ! This routines calculates the energy and the gradient of 
4
  ! a molecule, using VASP
5

    
6

    
7
  use Path_module, only : Nat, renum,Order,OrderInv, Coord, &
8
       ProgExe,Latr,FrozAtoms
9
  use Io_module
10

    
11
  !
12
  IMPLICIT NONE
13

    
14
  ! Energy (calculated if F300K=.F., else estimated)
15
  REAL(KREAL), INTENT (OUT) :: e
16
  ! Nb degree of freedom
17
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
18
  ! Gradient calculated at Geom geometry
19
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)
20

    
21
  ! ======================================================================
22

    
23
  character(LCHARS) :: LINE
24

    
25
  REAL(KREAL) :: Pi
26

    
27
  INTEGER(KINT) :: iat, i, n3at
28
  REAL(KREAL) :: RTmp1, RTmp2, RTmp3
29

    
30
  !
31
  CHARACTER(132), PARAMETER :: FileIn='POSCAR', FileOut='OUTCAR'
32

    
33
  CHARACTER(VLCHARS), SAVE :: RunCommand
34

    
35
!  REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:)
36

    
37
  LOGICAL :: Debug
38

    
39
  INTERFACE
40
     function valid(string) result (isValid)
41
       CHARACTER(*), intent(in) :: string
42
       logical                  :: isValid
43
     END function VALID
44
  END INTERFACE
45

    
46
  Pi=dacos(-1.0d0)
47
  debug=valid('EGRAD')
48
  if (debug) Call header('Entering Egrad_VASP')
49

    
50
  n3at=3*nat
51

    
52
!  ALLOCATE(GeomTmpC(n3at))
53

    
54

    
55
  gradcart=0.
56
  E=0.
57
  RunCommand=Trim(ProgExe) 
58

    
59
  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)
60

    
61
  ! we create the input file
62

    
63
!   write(*,*) "Coucou 2"
64

    
65
  OPEN(IOTMP,File=FileIn)
66

    
67
  ! we convert the coordinates into Cartesian coordinates
68
!     GeomTmpC=Reshape(GeomCart,SHAPE=(/3*Nat/))
69

    
70
   Call PrintGeomVasp(Vasp_Title,geomcart,Latr,FrozAtoms,IoTmp)
71

    
72
  CLOSE(IOTMP)
73

    
74

    
75
  IF (debug) THEN
76
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
77
     call system('cat ' // Trim(FileIn))
78
  END IF
79

    
80
  call system(RunCommand)
81

    
82
  if (debug) WRITE(*,*) 'DBG EGRAD_VASP, back from calculation'
83
!  if (debug) WRITE(*,*) "DBG EGRAD_VASP Order(I)",Order(1:Nat)
84

    
85
  ! Whatever the coordinate system, we read the cartesian gradient
86
  ! and we convert it to Zmat or other if needed
87

    
88
  OPEN(IOTMP,FILE=FileOut, STATUS='old')
89
  ! We first search for the forces
90
  READ(IOTMP,'(A)',END=999,ERR=999) LINE
91
! We search for the  last part of the OUTCAR file, after wavefunction is converged
92
  DO WHILE (INDEX(LINE,'EDIFF is reached')==0)
93
     READ(IOTMP,'(A)',END=999,ERR=999) LINE
94
  END DO
95
  READ(IOTMP,'(A)',END=999,ERR=999) LINE
96
  DO WHILE (INDEX(LINE,'TOTEN')==0)
97
     READ(IOTMP,'(A)',END=999,ERR=999) LINE
98
  END DO
99
  Line=Line(26:)
100
  READ(LINE,*) E
101
  if (debug) WRITE(*,'(1X,A,F20.6)') "DBG Egrad_VASP E=",E
102

    
103
! We search for the forces
104
  DO WHILE (INDEX(LINE,'TOTAL-FORCE')==0)
105
     READ(IOTMP,'(A)',END=999,ERR=999) LINE
106
  END DO
107
  READ(IOTMP,'(A)',END=999,ERR=999) LINE
108
  DO I=1,Nat
109
     Iat=I
110
     IF (renum) Iat=Order(I)
111
     READ(IOTMP,*,END=999,ERR=999)  Rtmp1,RTmp2,RTmp3,GradCart(3*Iat-2:3*Iat)
112
  END DO
113

    
114

    
115
  if (debug) THEN
116
     WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - original order"
117
     DO I=1,Nat
118
        Iat=Order(I)
119
        WRITE(*,'(1X,I5,3(1X,F20.6))') I,GradCart(3*Iat-2:3*Iat)
120
     END DO
121
     WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - internal order"
122
     DO I=1,Nat
123
        Iat=Order(I)
124
        WRITE(*,'(3(1X,I5),3(1X,F20.6))') I,Iat,OrderInv(I),GradCart(3*I-2:3*I)
125
     END DO
126
  END IF
127

    
128
  ! VASP gives the Forces in eV/Angstrom, so we convert it into the 
129
  ! gradient in ua/Angstrom 
130
    gradCart=-1.0d0*ev2Au*gradCart
131

    
132
  !  if (debug.AND.(COORD.EQ.'ZMAT')) THEN
133
  !     DO WHILE (INDEX(LINE,'Internal Coordinate Forces')==0)
134
  !        READ(IOTMP,'(A)') LINE
135
  !!        WRITE(*,*) 'Pas bon:' // TRIM(LINE)
136
  !     END DO
137
  !     WRITE(*,*) TRIM(LINE)
138
  !     DO I=1,Nat+2
139
  !        READ(IOTMP,'(A)') LINE
140
  !        WRITE(*,*) TRIM(LINE)
141
  !     END DO
142
  !  END IF
143

    
144
  CLOSE(IOTMP)
145

    
146
!  DeALLOCATE(GeomTmpC)
147

    
148

    
149
  if (debug) Call header('Egrad_VASP Over')
150

    
151
  RETURN
152

    
153
999 CONTINUE
154
  WRITE(*,*) 'We should not be here !!!!'
155
  STOP
156
  ! ======================================================================
157
end subroutine egrad_vasp
158

    
159

    
160

    
161