Statistiques
| Révision :

root / src / egrad_chamfre.f90 @ 8

Historique | Voir | Annoter | Télécharger (25,75 ko)

1
!*************############### Attention  ###############*************
2
!(a) Two parameters aa and br(3,3) must be provided in main program
3
!(b) At first, z=0 for the outmost surface 
4
!(c) The units of all parameters are not atomic units, but angstrom, eV etc.
5
!(d) Before using this routine, you must confirm that the energy is CONSISTENT with its derivative. 
6
!                        This check of consistency is VERY IMPORTANT for energy conservation in MD.
7
!!!
8
!!  xa, ya,za present the positions of atoms. v is the energy of the system. dvdr present the force along the x,y,z directions.
9
!!
10
subroutine egrad_chamfre(Nat,V,Geom,Grad)
11
! WAS xa,ya,za,v,dvdr)
12

    
13

    
14
  use Path_module, only : Lat_a,Lat_b,order,renum
15
  use Io_module
16

    
17

    
18
   IMPLICIT NONE
19

    
20
       INTEGER(KINT), INTENT(IN) :: Nat
21
       REAL(KREAL), INTENT(OUT) :: V,grad(Nat*3)
22
       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)
23

    
24
!     Bohr --> Angstr
25
!      real(KREAL), parameter :: BOHR   = 0.52917726D+00
26

    
27
!! PFL 2010 Nov. 28 ->
28
! for now I keep this old F77 scheme. This has to be replaced by ALLOCATABLE
29
! things
30
! <- PFL 2010 Nov 28.
31
  REAL(KREAL) :: xa(100),ya(100),za(100),dvdr(100,3)
32
  REAL(KREAL) :: rshort,rlong,zshort,zlong
33
  REAL(KREAL), SAVE :: rnc, rnn, rhh, rch, rnh, ann(100), anc(100), anh(100), ach(100), ahh(100)
34
  REAL(KREAL) :: dv(100,10,3),df(100,10,3)
35
  REAL(KREAL) :: px(9000),py(9000),pz(9000),tmp,temp,rtemp,yf
36
  INTEGER(KINT) :: ka(9000),kx(9000),ky(9000),itx,ity,nf(100),mf(100,500)
37
  REAL(KREAL) :: xi(3), xj(3), rij, er, eb
38
  REAL(KREAL) :: fcut, tf
39
  INTEGER(KINT) :: np,nq,nh,nc
40
  INTEGER(KINT) :: i, j, nn, m, n, mid, n1, im, jm, n2,n3,n4, iat
41
  LOGICAL, SAVE :: First=.TRUE.
42
  LOGICAL :: debug, FExist
43

    
44
  common/zdis/zshort,zlong
45
  common/rdistanc/rshort,rlong
46

    
47

    
48
  INTERFACE
49
     function valid(string) result (isValid)
50
       CHARACTER(*), intent(in) :: string
51
       logical                  :: isValid
52
     END function VALID
53
     
54
  END INTERFACE
55

    
56
!**********************************************
57

    
58

    
59
  debug=valid('egrad_chamfre')
60

    
61
 if (debug) Call Header('Entering Egrad_chamfre')
62

    
63
  v=0.0;dvdr=0.0
64
!******prepare the parameters of REBO potential
65
!rshort=3.5 
66
!rlong=3.9
67
zshort=3.5
68
zlong=4.5
69

    
70
rnn=3.2    !! to adjust
71
rnc=3.5    !! to adjust
72
rnh=3.0    !! to adjust
73
rch=3.2    !! to adjust
74
rhh=3.2    !! to adjust
75

    
76
!rph=rlong
77

    
78
 if (First) THEN
79
    First=.FALSE.
80
    INQUIRE(FILE="parameter.dat",EXIST=FExist)
81
    if (.NOT.Fexist) THEN
82
       WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!"
83
       WRITE(*,*) "!!   egrad_chamfre             !!!!"
84
       WRITE(*,*) "The file parameter.dat does not exists."
85
       WRITE(*,*) "!!   egrad_chamfre             !!!!"
86
       WRITE(*,*) "!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!"
87
       STOP
88
    END IF
89
open(1,file="parameter.dat")
90
do i=1,5   !Ni-Ni
91
read(1,*)ann(i)
92
enddo
93

    
94
do i=1,5   !Ni-C
95
read(1,*)anc(i)
96
enddo
97

    
98
do i=1,5   !Ni-H
99
read(1,*)anh(i)
100
enddo
101

    
102
do i=1,5   !H-H
103
read(1,*)ahh(i)
104
enddo
105
do i=1,5  !C-H
106
read(1,*)ach(i)
107
enddo
108
!********LONG INTERACTION **********
109
!do i=1,2
110
!read(1,*)a(i)
111
!enddo
112
close(1)
113
END IF
114
!************************************
115
np=45
116
nh=4
117
nc=1
118
nq=np+nh+nc
119
!************************************
120

    
121
 DO I=1,Nat
122
    Iat=Order(I)
123
    xa(I)=Geom(Iat,1)
124
    ya(I)=Geom(Iat,2)
125
    za(I)=Geom(Iat,3)
126
 END DO
127

    
128
px=0.0;py=0.0;pz=0.0
129
nn=0
130
 do m=-1,1
131
    do n=-1,1
132
       do i=1,nq
133
          if(m==0.and.n==0.and.i==1)mid=nn
134
            nn=nn+1
135
            px(nn)=xa(i)+m*lat_a(1)+n*lat_b(1)
136
            py(nn)=ya(i)+m*lat_a(2)+n*lat_b(2)
137
            pz(nn)=za(i)
138
            ka(nn)=i;kx(nn)=m;ky(nn)=n
139
       enddo
140
    enddo
141
 enddo
142
!***********************xyz file required for MS input
143
 open(123,file="111.xyz",ACCESS='APPEND')
144
    write(123,*)nn
145
    write(123,*)
146
   do i=1,nn
147
       if(ka(i)<(np+1))then
148
            write(123,10)'Ni',px(i),py(i),pz(i)
149
       elseif (ka(i)<np+5) then
150
            write(123,10)'H',px(i),py(i),pz(i)
151
       else
152
            write(123,10)'C',px(i),py(i),pz(i)
153
  10 format(a6,3f15.8)
154
       endif
155
   enddo
156
 close(123)
157
!******************************************
158
 do i=1,nq
159
    tmp=px(i);px(i)=px(i+mid);px(i+mid)=tmp
160
    tmp=py(i);py(i)=py(i+mid);py(i+mid)=tmp
161
    tmp=pz(i);pz(i)=pz(i+mid);pz(i+mid)=tmp
162
    itx=kx(i);kx(i)=kx(i+mid);kx(i+mid)=itx
163
    ity=ky(i);ky(i)=ky(i+mid);ky(i+mid)=ity
164
 enddo
165
 mid=0
166
!******************************************
167
  nf=0;mf=0
168
    do i=1,nq
169
        nf(i)=0
170
           do j=1,nn
171
             temp=sqrt((px(i)-px(j))**2+(py(i)-py(j))**2+(pz(i)-pz(j))**2)
172
!
173
!rtemp=rlong
174
    if (ka(i+mid)<np+1.and.ka(j+mid)<np+1)then     ! Ni-Ni
175
       
176
        rtemp=rnn
177
   endif
178
    if (ka(i+mid)<np+1.and.ka(j+mid)>np.and.ka(j+mid)<np+5)then  ! Ni-H
179
          
180
           rtemp=rnh
181
     endif
182
    if (ka(i+mid)<np+1.and.ka(j+mid)>np+4) then   ! Ni-C
183
           
184
           rtemp=rnc
185
     endif
186
      if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)>np+4)  then ! H-C
187
           
188
            rtemp=rch
189
     endif
190
    if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)>np.and.ka(j+mid)<np+5)then  ! H-H
191
         
192
          rtemp=rhh
193
   endif
194
    if (ka(i+mid)>np.and.ka(i+mid)<np+5.and.ka(j+mid)<np+1) then  ! H-Ni
195
            
196
             rtemp=rnh
197
       endif
198
    if (ka(i+mid)>np+4.and.ka(j+mid)<np+1) then !  C-Ni
199
            
200
             rtemp=rnc
201
     endif
202
    if (ka(i+mid)>np+4.and.ka(j+mid)>np.and.ka(j+mid)<np+5) then ! C-H
203
             
204
              rtemp=rch
205
     endif
206

    
207
      if(temp>0.1.and.temp<rtemp)then
208
         nf(i)=nf(i)+1
209
         mf(i,nf(i))=j
210
      endif
211
!
212
    enddo
213
  !write(*,*)nf(i)
214
 enddo
215
!**************************************
216
!stop
217

    
218
!*************************************Starting to calculate the sum
219
  yf=0.0;dvdr=0.0
220
   do i=1,nq
221
      er=0.0;eb=0.0
222
      dv=0.0;df=0.0
223
      do j=1,nf(i)
224

    
225
             xi(1)=px(i);      xi(2)=py(i);      xi(3)=pz(i)
226
             xj(1)=px(mf(i,j));xj(2)=py(mf(i,j));xj(3)=pz(mf(i,j))
227
             rij=sqrt((xj(1)-xi(1))**2.0+(xj(2)-xi(2))**2.0+(xj(3)-xi(3))**2.0)
228
             n1=ka(i)/(np+1);n2=ka(mf(i,j))/(np+1)    !!!  Ni 46  H  50  C
229
             n3=ka(i)/(np+5);n4=ka(mf(i,j))/(np+5)
230
  
231
 !**********************************Ni-Ni interaction
232
 ! **********
233
         if(n1==0.and.n3==0.and.n2==0.and.n4==0)then     !Ni-Ni
234
           call fctNiNi(rnn,rij,fcut)
235
           call dfcutNiNi(rnn,rij,tf)
236
            eb=eb+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
237
            er=er+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut
238

    
239
!Calculation of FORCE
240
do im=1,3
241
   dv(i,1,im)=dv(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*fcut &
242
        *(-ann(2)/ann(3))*(xi(im)-xj(im))/rij
243
   dv(i,2,im)=dv(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
244
   dv(i,3,im)=dv(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut &
245
        *(-2.0*ann(1)/ann(3))*(xi(im)-xj(im))/rij
246
   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2) &
247
     *(rij/ann(3)-1.0))*fcut*(-ann(2)/ann(3))*(xj(im)-xi(im))/rij
248
   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1) &
249
     *(rij/ann(3)-1.0))*fcut
250
   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1) &
251
        *(rij/ann(3)-1.0))*fcut*(-2.0*ann(1)/ann(3))*(xj(im)-xi(im))/rij
252

    
253
   df(i,1,im)=df(i,1,im)+ann(4)*exp(-1.0*ann(2)*(rij/ann(3)-1.0))*tf &
254
        *(xi(im)-xj(im))/rij
255
   df(i,2,im)=df(i,2,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*fcut
256
   df(i,3,im)=df(i,3,im)+ann(5)*exp(-2.0*ann(1)*(rij/ann(3)-1.0))*tf &
257
        *(xi(im)-xj(im))/rij
258
   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ann(4)*exp(-1.0*ann(2) &
259
     *(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij
260
   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ann(5)*exp(-2.0*ann(1) &
261
        *(rij/ann(3)-1.0))*fcut
262
   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ann(5)*exp(-2.0*ann(1) &
263
        *(rij/ann(3)-1.0))*tf*(xj(im)-xi(im))/rij
264
enddo
265
!
266
endif
267

    
268
 !**********************************Ni-H interaction
269
! **********
270
         if(n1==0.and.n3==0.and.n2==1.and.n4==0)then     !Ni-H
271
           call fctNiH(rnh,rij,fcut)
272
           call dfcutNiH(rnh,rij,tf)
273
             eb=eb+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
274
             er=er+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut
275

    
276
!Calculation of FORCE
277
do im=1,3
278
   dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut &
279
        *(-anh(2)/anh(3))*(xi(im)-xj(im))/rij
280
   dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
281
   dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut &
282
        *(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij
283
   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)* &
284
        exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3)) &
285
        *(xj(im)-xi(im))/rij
286
   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)* &
287
        exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
288
   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)* &
289
        exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3)) &
290
        *(xj(im)-xi(im))/rij
291

    
292
   df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf &
293
        *(xi(im)-xj(im))/rij
294
   df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
295
   df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf &
296
        *(xi(im)-xj(im))/rij
297
   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2) &
298
        *(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
299
   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0*anh(1) &
300
        *(rij/anh(3)-1.0))*fcut
301
   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0*anh(1) &
302
        *(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
303
enddo
304
!
305
!
306
             endif
307

    
308
!**********************************Ni-C interaction
309
! **********
310
         if(n1==0.and.n3==0.and.n2==1.and.n4==1)then     !Ni-C
311
           call fctNiC(rnc,rij,fcut)
312
           call dfcutNiC(rnc,rij,tf)
313
            eb=eb+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
314
            er=er+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut
315

    
316
!Calculation of FORCE
317
do im=1,3
318
   dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut &
319
        *(-anc(2)/anc(3))*(xi(im)-xj(im))/rij
320
   dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
321
   dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut &
322
        *(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij
323
   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
324
        anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij
325
   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
326
        anc(1)*(rij/anc(3)-1.0))*fcut
327
   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
328
        anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij
329
   
330
   df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf &
331
        *(xi(im)-xj(im))/rij
332
   df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
333
   df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf &
334
        *(xi(im)-xj(im))/rij
335
   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
336
        anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
337
   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
338
        anc(1)*(rij/anc(3)-1.0))*fcut
339
   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
340
        anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
341
enddo
342
!
343
               endif
344
!!**********************************H-Ni interaction
345
! **********
346
         if(n1==1.and.n3==0.and.n2==0.and.n4==0)then        ! H-Ni
347
           call fctNiH(rnh,rij,fcut)
348
           call dfcutNiH(rnh,rij,tf)
349
             eb=eb+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
350
             er=er+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut
351

    
352
!Calculation of FORCE
353
do im=1,3
354
   dv(i,1,im)=dv(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*fcut &
355
        *(-anh(2)/anh(3))*(xi(im)-xj(im))/rij
356
   dv(i,2,im)=dv(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
357
   dv(i,3,im)=dv(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut &
358
        *(-2.0*anh(1)/anh(3))*(xi(im)-xj(im))/rij
359
   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0*anh(2) &
360
     *(rij/anh(3)-1.0))*fcut*(-anh(2)/anh(3))*(xj(im)-xi(im))/rij
361
   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0* & 
362
        anh(1)*(rij/anh(3)-1.0))*fcut
363
   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0* &
364
        anh(1)*(rij/anh(3)-1.0))*fcut*(-2.0*anh(1)/anh(3))*(xj(im)-xi(im))/rij
365

    
366
   df(i,1,im)=df(i,1,im)+anh(4)*exp(-1.0*anh(2)*(rij/anh(3)-1.0))*tf &
367
        *(xi(im)-xj(im))/rij
368
   df(i,2,im)=df(i,2,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*fcut
369
   df(i,3,im)=df(i,3,im)+anh(5)*exp(-2.0*anh(1)*(rij/anh(3)-1.0))*tf &
370
        *(xi(im)-xj(im))/rij
371
   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anh(4)*exp(-1.0* &
372
        anh(2)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
373
   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anh(5)*exp(-2.0* &
374
        anh(1)*(rij/anh(3)-1.0))*fcut
375
   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anh(5)*exp(-2.0* &
376
        anh(1)*(rij/anh(3)-1.0))*tf*(xj(im)-xi(im))/rij
377
enddo
378
!
379
              endif
380
!!**********************************H-H interaction
381
! **********
382
        if(n1==1.and.n3==0.and.n2==1.and.n4==0)then   !H-H
383
           call fctHH(rhh,rij,fcut)
384
           call dfcutHH(rhh,rij,tf)
385
              eb=eb+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
386
              er=er+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut
387

    
388
!Calculation of FORCE
389
do im=1,3
390
   dv(i,1,im)=dv(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*fcut &
391
        *(-ahh(2)/ahh(3))*(xi(im)-xj(im))/rij
392
   dv(i,2,im)=dv(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
393
   dv(i,3,im)=dv(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut &
394
       *(-2.0*ahh(1)/ahh(3))*(xi(im)-xj(im))/rij
395
   dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0* &
396
       ahh(2)*(rij/ahh(3)-1.0))*fcut*(-ahh(2)/ahh(3))*(xj(im)-xi(im))/rij
397
   dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0* &
398
       ahh(1)*(rij/ahh(3)-1.0))*fcut
399
   dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0* &
400
       ahh(1)*(rij/ahh(3)-1.0))*fcut*(-2.0*ahh(1)/ahh(3))*(xj(im)-xi(im))/rij
401
   
402
   df(i,1,im)=df(i,1,im)+ahh(4)*exp(-1.0*ahh(2)*(rij/ahh(3)-1.0))*tf &
403
        *(xi(im)-xj(im))/rij
404
   df(i,2,im)=df(i,2,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*fcut
405
   df(i,3,im)=df(i,3,im)+ahh(5)*exp(-2.0*ahh(1)*(rij/ahh(3)-1.0))*tf &
406
        *(xi(im)-xj(im))/rij
407
   df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ahh(4)*exp(-1.0* &
408
        ahh(2)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij
409
   df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ahh(5)*exp(-2.0* &
410
        ahh(1)*(rij/ahh(3)-1.0))*fcut
411
   df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ahh(5)*exp(-2.0* &
412
        ahh(1)*(rij/ahh(3)-1.0))*tf*(xj(im)-xi(im))/rij
413
enddo
414
!
415
             endif
416

    
417

    
418
!!**********************************H-C interaction
419
! **********
420
        if(n1==1.and.n3==0.and.n2==1.and.n4==1)then   !H-C
421
            call fctCH(rch,rij,fcut)
422
            call dfcutCH(rch,rij,tf)
423
             eb=eb+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
424
             er=er+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut
425

    
426
!Calculation of FORCE
427
do im=1,3
428
dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut &
429
        *(-ach(2)/ach(3))*(xi(im)-xj(im))/rij
430
dv(i,2,im)=dv(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
431
dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut &
432
        *(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij
433
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
434
        ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij
435
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
436
        ach(1)*(rij/ach(3)-1.0))*fcut
437
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
438
        ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij
439

    
440
df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf &
441
        *(xi(im)-xj(im))/rij
442
df(i,2,im)=df(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
443
df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf &
444
        *(xi(im)-xj(im))/rij
445
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
446
        ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
447
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
448
        ach(1)*(rij/ach(3)-1.0))*fcut
449
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
450
        ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
451
enddo
452
!
453
                endif
454
!!**********************************C-Ni interaction
455
! **********
456
        if(n1==1.and.n3==1.and.n2==0.and.n4==0)then   !C-Ni
457
            call fctNiC(rnc,rij,fcut)
458
            call dfcutNiC(rnc,rij,tf)
459
             eb=eb+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
460
             er=er+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut
461

    
462
!Calculation of FORCE
463
do im=1,3
464
dv(i,1,im)=dv(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*fcut &
465
        *(-anc(2)/anc(3))*(xi(im)-xj(im))/rij
466
dv(i,2,im)=dv(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
467
dv(i,3,im)=dv(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut &
468
        *(-2.0*anc(1)/anc(3))*(xi(im)-xj(im))/rij
469
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
470
        anc(2)*(rij/anc(3)-1.0))*fcut*(-anc(2)/anc(3))*(xj(im)-xi(im))/rij
471
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
472
        anc(1)*(rij/anc(3)-1.0))*fcut
473
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
474
        anc(1)*(rij/anc(3)-1.0))*fcut*(-2.0*anc(1)/anc(3))*(xj(im)-xi(im))/rij
475

    
476
df(i,1,im)=df(i,1,im)+anc(4)*exp(-1.0*anc(2)*(rij/anc(3)-1.0))*tf &
477
        *(xi(im)-xj(im))/rij
478
df(i,2,im)=df(i,2,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*fcut
479
df(i,3,im)=df(i,3,im)+anc(5)*exp(-2.0*anc(1)*(rij/anc(3)-1.0))*tf &
480
        *(xi(im)-xj(im))/rij
481
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+anc(4)*exp(-1.0* &
482
        anc(2)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
483
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+anc(5)*exp(-2.0* &
484
        anc(1)*(rij/anc(3)-1.0))*fcut
485
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+anc(5)*exp(-2.0* &
486
        anc(1)*(rij/anc(3)-1.0))*tf*(xj(im)-xi(im))/rij
487
enddo
488
!
489

    
490
               endif
491
!!**********************************C-H interaction
492
! **********
493
       if(n1==1.and.n3==1.and.n2==1.and.n4==0)then   !C-H
494
           call fctCH(rch,rij,fcut)
495
           call dfcutCH(rch,rij,tf)
496
             eb=eb+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
497
             er=er+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut
498

    
499

    
500
   !Calculation of FORCE
501
do im=1,3
502
dv(i,1,im)=dv(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*fcut &
503
       *(-ach(2)/ach(3))*(xi(im)-xj(im))/rij
504
dv(i,2,im)=dv(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
505
dv(i,3,im)=dv(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut &
506
       *(-2.0*ach(1)/ach(3))*(xi(im)-xj(im))/rij
507
dv(ka(mf(i+mid,j)),1,im)=dv(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
508
       ach(2)*(rij/ach(3)-1.0))*fcut*(-ach(2)/ach(3))*(xj(im)-xi(im))/rij
509
dv(ka(mf(i+mid,j)),2,im)=dv(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
510
       ach(1)*(rij/ach(3)-1.0))*fcut
511
dv(ka(mf(i+mid,j)),3,im)=dv(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
512
       ach(1)*(rij/ach(3)-1.0))*fcut*(-2.0*ach(1)/ach(3))*(xj(im)-xi(im))/rij
513

    
514
df(i,1,im)=df(i,1,im)+ach(4)*exp(-1.0*ach(2)*(rij/ach(3)-1.0))*tf &
515
       *(xi(im)-xj(im))/rij
516
df(i,2,im)=df(i,2,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*fcut
517
df(i,3,im)=df(i,3,im)+ach(5)*exp(-2.0*ach(1)*(rij/ach(3)-1.0))*tf &
518
       *(xi(im)-xj(im))/rij
519
df(ka(mf(i+mid,j)),1,im)=df(ka(mf(i+mid,j)),1,im)+ach(4)*exp(-1.0* &
520
       ach(2)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
521
df(ka(mf(i+mid,j)),2,im)=df(ka(mf(i+mid,j)),2,im)+ach(5)*exp(-2.0* &
522
       ach(1)*(rij/ach(3)-1.0))*fcut
523
df(ka(mf(i+mid,j)),3,im)=df(ka(mf(i+mid,j)),3,im)+ach(5)*exp(-2.0* &
524
       ach(1)*(rij/ach(3)-1.0))*tf*(xj(im)-xi(im))/rij
525
enddo
526
!
527

    
528
        endif
529
  !!***************************ONE end****************      
530

    
531
enddo
532
!!******** CALCULATION ENERGY
533
!
534
        yf=yf+er-sqrt(eb)
535
!!*********
536
       do im=1,3
537
            dvdr(i,im)=dvdr(i,im)+dv(i,1,im)-0.5/(dv(i,2,im))**0.5*dv(i,3,im)
538
            dvdr(i,im)=dvdr(i,im)+df(i,1,im)-0.5/(df(i,2,im))**0.5*df(i,3,im)
539
         do j=1,nf(i+mid)
540
             jm=ka(mf(i+mid,j))
541
             dvdr(jm,im)=dvdr(jm,im)+dv(jm,1,im)-0.5/(dv(i,2,im))**0.5*dv(jm,3,im)
542
             dvdr(jm,im)=dvdr(jm,im)+df(jm,1,im)-0.5/(df(i,2,im))**0.5*df(jm,3,im)
543
         enddo
544
       enddo
545
!
546
enddo
547
!!*****************************TWO ENDS**************************
548
!***************************************Long-ranged interaction
549
!do i=nq-4,nq
550
!     zi=pz(i+mid) !-0.3*aa*br(3,3)
551
!     call fct2(zshort,zlong,zi,fcut2)
552
 !    call dfcut2(zshort,zlong,zi,tf)
553
 !    yf=yf+fcut2*(a(6)-a(7)/zi**2.0)
554

    
555
 !  dvdr(i,3)=dvdr(i,3)+fcut2*(a(7)*2.0/zi**3.0)+tf*(a(6)-a(7)/zi**2.0)
556

    
557
!enddo
558
!*******************************************************************
559
  v=yf 
560
   do im=1,3
561
      do jm=1,nq
562
          dvdr(jm,im)=dvdr(jm,im)  
563
      enddo
564
  enddo
565

    
566
!PFL: We convert dvdr into Grad
567
      do jm=1,nq
568
          Iat=Jm
569
          if (renum) Iat=order(jm)
570
          Grad(3*Iat-2:3*Iat)=dvdr(jm,1:3)
571
      enddo
572
  
573
 if (debug) WRITE(*,*) "EGRAD_CHAMFRE, E=",V
574
  
575
 if (debug) call Header('Egrad_chamfre OVER')
576

    
577
!********************************************ENDENDENDENDENDENDENDEND
578
return
579
end subroutine egrad_chamfre
580

    
581
!******************************************
582
!**********PARAMETERS**********************
583
subroutine fctNiNi(r2,r,fcut)
584

    
585
   IMPLICIT NONE
586
!     integer, parameter :: KINT = kind(1)
587
     integer, parameter :: KREAL = kind(1.0d0)
588

    
589
 REAL(KREAL) :: r2,r,fcut,r1,pi
590
pi=4.0*atan(1.0)
591
r1=r2-0.5
592

    
593
if(r<=r1)then
594
fcut=1.0
595
endif
596
if(r>r1.and.r<=r2)then
597
fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
598
endif
599
if(r>r2)then
600
fcut=0.0
601
endif
602
return
603
end subroutine fctNiNi
604
!!************
605
subroutine dfcutNiNi(r2,r,df)
606
   IMPLICIT NONE
607
!     integer, parameter :: KINT = kind(1)
608
     integer, parameter :: KREAL = kind(1.0d0)
609
  REAL(KREAL) :: r, r2,df,pi,r1
610
pi=4.0*atan(1.0)
611
r1=r2-0.5
612

    
613
if(r<=r1)then
614
df=0.0
615
endif
616
if(r>r1.and.r<=r2)then
617
df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
618
endif
619
if(r>r2)then
620
df=0.0
621
endif
622
return
623
end subroutine dfcutNiNi
624

    
625

    
626
!!************
627
subroutine fctNiH(r2,r,fcut)
628
   IMPLICIT NONE
629
!     integer, parameter :: KINT = kind(1)
630
     integer, parameter :: KREAL = kind(1.0d0)
631

    
632
 REAL(KREAL) ::r2,r,fcut
633
REAL(KREAL) ::rshort,rlong,r1,pi
634

    
635
 rshort=2.7
636
 rlong=3.0
637
 pi=4.0*atan(1.0)
638
 r1=rshort
639
 r2=rlong
640
if(r<=r1)then
641
 fcut=1.0
642
endif
643
if(r>r1.and.r<=r2)then
644
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
645
endif
646
if(r>r2)then
647
 fcut=0.0
648
endif
649
return
650
end subroutine fctNiH
651
!!**********
652

    
653
subroutine dfcutNiH(r2,r,df)
654

    
655

    
656
   IMPLICIT NONE
657
!     integer, parameter :: KINT = kind(1)
658
     integer, parameter :: KREAL = kind(1.0d0)
659

    
660

    
661
REAL(KREAL) ::r2,r,df
662
REAL(KREAL) ::rshort,rlong,r1,pi
663

    
664
 rshort=2.7
665
 rlong=3.0
666
 pi=4.0*atan(1.0)
667
 r1=rshort
668
 r2=rlong
669
if(r<=r1)then
670
 df=0.0
671
endif
672
if(r>r1.and.r<=r2)then
673
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
674
endif
675
if(r>r2)then
676
 df=0.0
677
endif
678
return
679
end subroutine dfcutNiH
680

    
681
!!************
682
subroutine fctNiC(r2,r,fcut)
683

    
684
   IMPLICIT NONE
685
!     integer, parameter :: KINT = kind(1)
686
     integer, parameter :: KREAL = kind(1.0d0)
687

    
688

    
689
REAL(KREAL) ::r2,r,fcut
690
REAL(KREAL) ::rshort,rlong,r1,pi
691

    
692
 rshort=3.2
693
 rlong=3.5
694
 pi=4.0*atan(1.0)
695
 r1=rshort
696
 r2=rlong
697
if(r<=r1)then
698
 fcut=1.0
699
endif
700
if(r>r1.and.r<=r2)then
701
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
702
endif
703
if(r>r2)then
704
 fcut=0.0
705
endif
706
return
707
end subroutine fctNiC
708

    
709
!!***********
710
subroutine dfcutNiC(r2,r,df)
711

    
712
   IMPLICIT NONE
713
!     integer, parameter :: KINT = kind(1)
714
     integer, parameter :: KREAL = kind(1.0d0)
715

    
716
     
717
REAL(KREAL) ::r2,r,df
718
REAL(KREAL) ::rshort,rlong,r1,pi
719

    
720
 rshort=3.2
721
 rlong=3.5
722
 pi=4.0*atan(1.0)
723
 r1=rshort
724
 r2=rlong
725
if(r<=r1)then
726
 df=0.0
727
endif
728
if(r>r1.and.r<=r2)then
729
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
730
endif
731
if(r>r2)then
732
 df=0.0
733
endif
734
return
735
end subroutine dfcutNiC
736

    
737
!!************
738
subroutine fctCH(r2,r,fcut)
739

    
740

    
741
   IMPLICIT NONE
742
!     integer, parameter :: KINT = kind(1)
743
     integer, parameter :: KREAL = kind(1.0d0)
744

    
745

    
746
REAL(KREAL) ::r2,r,fcut
747
REAL(KREAL) ::rshort,rlong,r1,pi
748

    
749
 rshort=2.1
750
 rlong=3.2
751
 pi=4.0*atan(1.0)
752
 r1=rshort
753
 r2=rlong
754
if(r<=r1)then
755
 fcut=1.0
756
endif
757
if(r>r1.and.r<=r2)then
758
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
759
endif
760
if(r>r2)then
761
 fcut=0.0
762
endif
763
return
764
end subroutine fctCH
765
!!*****************
766

    
767
subroutine dfcutCH(r2,r,df)
768

    
769

    
770
   IMPLICIT NONE
771
!     integer, parameter :: KINT = kind(1)
772
     integer, parameter :: KREAL = kind(1.0d0)
773

    
774
     
775
REAL(KREAL) ::r2,r,df
776
REAL(KREAL) ::rshort,rlong,r1,pi
777

    
778
 rshort=2.1
779
 rlong=3.2
780
 pi=4.0*atan(1.0)
781
 r1=rshort
782
 r2=rlong
783
if(r<=r1)then
784
 df=0.0
785
endif
786
if(r>r1.and.r<=r2)then
787
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
788
endif
789
if(r>r2)then
790
 df=0.0
791
endif
792
return
793
end subroutine dfcutCH
794
!!************
795
subroutine fctHH(r2,r,fcut)
796

    
797
   IMPLICIT NONE
798
!     integer, parameter :: KINT = kind(1)
799
     integer, parameter :: KREAL = kind(1.0d0)
800

    
801

    
802
REAL(KREAL) ::r2,r,fcut,r1,pi
803
pi=4.0*atan(1.0)
804
r1=r2-1.1
805
if(r<=r1)then
806
 fcut=1.0
807
endif
808
if(r>r1.and.r<=r2)then
809
 fcut=(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
810
endif
811
if(r>r2)then
812
 fcut=0.0
813
endif
814
return
815
end subroutine fctHH
816
!!************
817
subroutine dfcutHH(r2,r,df)
818

    
819

    
820
   IMPLICIT NONE
821
!     integer, parameter :: KINT = kind(1)
822
     integer, parameter :: KREAL = kind(1.0d0)
823

    
824

    
825
REAL(KREAL) ::r2,r,df,r1,pi
826
pi=4.0*atan(1.0)
827
r1=r2-1.1
828
if(r<=r1)then
829
 df=0.0
830
endif
831
if(r>r1.and.r<=r2)then
832
 df=-pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
833
endif
834
if(r>r2)then
835
 df=0.0
836
endif
837
return
838
end subroutine dfcutHH
839

    
840
!!*********************LONG INTERACTION
841

    
842
subroutine dfcut2(zshort,zlong,r,df)
843

    
844

    
845
   IMPLICIT NONE
846
!     integer, parameter :: KINT = kind(1)
847
     integer, parameter :: KREAL = kind(1.0d0)
848

    
849

    
850
REAL(KREAL) ::zshort,zlong,r,df,pi, r1, r2
851

    
852
pi=4.0*atan(1.0)
853

    
854
r1=zshort
855
r2=zlong
856

    
857
if(r<=r1)then
858
df=0.0
859
endif
860
if(r>r1.and.r<=r2)then
861
df=pi/2.0/(r2-r1)*sin(pi*(r-r1)/(r2-r1))
862
endif
863
if(r>r2)then
864
df=0.0
865
endif
866

    
867
return
868
end subroutine dfcut2
869

    
870
subroutine fct2(zshort,zlong,r,fcut)
871

    
872
   IMPLICIT NONE
873
!     integer, parameter :: KINT = kind(1)
874
     integer, parameter :: KREAL = kind(1.0d0)
875

    
876

    
877
REAL(KREAL) ::zshort,zlong,r,fcut,pi,r1,r2
878

    
879
pi=4.0*atan(1.0)
880

    
881
r1=zshort
882
r2=zlong
883

    
884
if(r<=r1)then
885
fcut=0.0
886
endif
887
if(r>r1.and.r<=r2)then
888
fcut=1.0-(1.0+cos(pi*(r-r1)/(r2-r1)))/2.0
889
endif
890
if(r>r2)then
891
fcut=1.0
892
endif
893

    
894
return
895
end subroutine fct2
896