c write(*,*) 'case1:',tmp_n
- t1dx=t1dx+0.0d0
- t1dy=t1dy+0.0d0
- t1dz=t1dz+0.0d0
+cc t1dx=t1dx+0.0d0
+cc t1dy=t1dy+0.0d0
+cc t1dz=t1dz+0.0d0
t2dx(j)=0.0d0
t2dy(j)=0.0d0
t2dz(j)=0.0d0
tmp_n=tmp_n+1.0d0
c write(*,*) 'case4:',tmp_n
- t1dx=t1dx+0.0d0
- t1dy=t1dy+0.0d0
- t1dz=t1dz+0.0d0
+cc t1dx=t1dx+0.0d0
+cc t1dy=t1dy+0.0d0
+cc t1dz=t1dz+0.0d0
t2dx(j)=0.0d0
t2dy(j)=0.0d0
t2dz(j)=0.0d0
double precision enebet
enebet=0.0d0
- bx=0.0d0;by=0.0d0;bz=0.0d0
- shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
+c bx=0.0d0;by=0.0d0;bz=0.0d0
+c shetfx=0.0d0;shetfy=0.0d0;shetfz=0.0d0
gdfab=0.0d0
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc******************************************************************************
c real*8 dfaexp(15001)
real*8 uum, uup
real*8 vbeta,vbetp,vbetm,y,y1,y2,yshe1,yshe2,yy1,yy2
+ real*8 shetfx(maxca),shetfy(maxca),shetfz(maxca)
+ common /shetf/ shetfx,shetfy,shetfz
+
common /sheca/ bx,by,bz
common /phys1/ inb,nmax,iselect
common /kyori2/ dis
common /sinsin/ sth
real*8 r_pair_mat(maxca,maxca)
+ real*8 e_gcont,fprim_gcont,de_gcont
ci integer istrand(maxca,maxca)
ci integer istrand_p(maxca,maxca),istrand_m(maxca,maxca)
ci common /shetest/ istrand,istrand_p,istrand_m
do i=1,inb-7
do j=i+4,inb-3
+
+ if (dis(i,j).lt.dfa_cutoff) then
+ call gcont(dis(i,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
+
ip=i+1
ipp=i+2
jp=j+1
vbet(i,j)=uup+uum
vbetp=vbetp+uup
vbetm=vbetm+uum
- vbeta=vbeta+vbet(i,j)
+ vbeta=vbeta+vbet(i,j)*e_gcont
+
+
+ if (dis(i,j) .ge. dfa_cutoff-2*dfa_cutoff_delta) then
+c gradient correction from gcont
+ de_gcont=vbet(i,j)*fprim_gcont/dis(i,j)
+ shetfx(i)=shetfx(i) + de_gcont*rx(i,j)
+ shetfy(i)=shetfy(i) + de_gcont*ry(i,j)
+ shetfz(i)=shetfz(i) + de_gcont*rz(i,j)
+
+ shetfx(j)=shetfx(j) - de_gcont*rx(i,j)
+ shetfy(j)=shetfy(j) - de_gcont*ry(i,j)
+ shetfz(j)=shetfz(j) - de_gcont*rz(i,j)
+
+c energy correction from gcont
+ vbet(i,j)=vbet(i,j)*e_gcont
+ vbetap(i,j)=vbetap(i,j)*e_gcont
+ vbetap1(i,j)=vbetap1(i,j)*e_gcont
+ vbetap2(i,j)=vbetap2(i,j)*e_gcont
+ vbetam(i,j)=vbetam(i,j)*e_gcont
+ vbetam1(i,j)=vbetam1(i,j)*e_gcont
+ vbetam2(i,j)=vbetam2(i,j)*e_gcont
+ endif
+
ci elseif(istrand(i,j).eq.0)then
ci vbet(i,j)=0
c write(*,*) 'vbetm:',vbetm
c write(*,*) 'vbet(i,j):',vbet(i,j)
c stop
-
+
+ else
+ vbetap(i,j)=0
+ vbetap1(i,j)=0
+ vbetap2(i,j)=0
+ vbetam(i,j)=0
+ vbetam1(i,j)=0
+ vbetam2(i,j)=0
+ vbet(i,j)=0
+ endif
enddo
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbet(maxca,maxca)
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 c3,v3,cc2,cc3,dmm3,dmm3__,c4,v4,c7,v7,cc7,c8,v8,cc8
real*8 c9,v9,cc9,dmm9,dmm9__,c10,v10,dmm2,dmm1__,dmm2_1,dmm2_2
real*8 dmm7,dmm8,dmm7__,dmm8_1,dmm8_2
+ real*8 e_gcont,fprim_gcont
C--------------------------------------------------------------------------------
do i=4,inb-4
im3=i-3
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 y66x,y66y,y66z,yy6,yyy4,yyy5a,yyy5b
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b
+ real*8 e_gcont,fprim_gcont
c********************************************************************************
do i=3,inb-5
imm=i-2
im=i-1
do j=i+2,inb-3
+
+ if (dis(imm,j).lt.dfa_cutoff) then
+ call gcont(dis(imm,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
jp=j+1
jpp=j+2
! $ -sy1*vbetam1(imm,j)-sy2*vbetam2(imm,j)
! shefz(i,5)=shefz(i,5)
! $ -sz1*vbetam1(imm,j)-sz2*vbetam2(imm,j)
-
+
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy8,yyy9a,yyy9b,yyy4
real*8 yyy3a,yyy3b,y66z,yy6,yyy5a,yyy5b
+ real*8 e_gcont,fprim_gcont
C********************************************************************************
do i=2,inb-6
ip=i+1
im=i-1
do j=i+3,inb-3
+
+ if (dis(im,j).lt.dfa_cutoff) then
+ call gcont(dis(im,j),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
jp=j+1
jpp=j+2
! $ -sy1*vbetam1(im,j)-sy2*vbetam2(im,j)
! shefz(i,6)=shefz(i,6)
! $ -sz1*vbetam1(im,j)-sz2*vbetam2(im,j)
-
+
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy4,yyy5a,yyy5b,yy6
real*8 yyy9a,yyy9b,y5z,y66z,y9z,yyy8
+ real*8 e_gcont,fprim_gcont
C********************************************************************************
do j=7,inb-1
jm=j-1
jmm=j-2
do i=1,j-6
+
+ if (dis(i,jmm).lt.dfa_cutoff) then
+ call gcont(dis(i,jmm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
ip=i+1
ipp=i+2
! shefz(j,11)=shefz(j,11)
! $ -sz1*vbetam1(i,jmm)-sz2*vbetam2(i,jmm)
+ endif
ci endif
enddo
implicit none
integer maxca
parameter(maxca=800)
+ real*8 dfa_cutoff,dfa_cutoff_delta
+ parameter(dfa_cutoff=15.5d0)
+ parameter(dfa_cutoff_delta=0.5d0)
cc**********************************************************************
real*8 vbetap(maxca,maxca),vbetam(maxca,maxca)
real*8 vbetap1(maxca,maxca),vbetam1(maxca,maxca)
real*8 sx,sy,sz,sx1,sy1,sz1,sx2,sy2,sz2,y6x,y6y,y6z,y66x,y66y,y66z
real*8 yy88,yyy8a,yyy8b,yy8,y8x,y8y,y8z,yy99,yyy9,yy9,y9x,y9y,y9z
real*8 yy1010,yyy10,yy10,y10x,y10y,y10z,yyy10a,yyy10b,yy6,yyy8
+ real*8 e_gcont,fprim_gcont
!c*************************************************************************c
do j=6,inb-2
jp=j+1
jm=j-1
do i=1,j-5
+
+ if (dis(i,jm).lt.dfa_cutoff) then
+ call gcont(dis(i,jm),dfa_cutoff-dfa_cutoff_delta,1.0D0,
+ & dfa_cutoff_delta,e_gcont,fprim_gcont)
+
ip=i+1
ipp=i+2
shefz(j,12)=shefz(j,12)-sz*vbetam(i,jm)
$ -sz1*vbetam1(i,jm)-sz2*vbetam2(i,jm)
+ endif
+
ci endif
ENDDO