include 'COMMON.IOUNITS'
include 'COMMON.CHAIN'
include 'COMMON.DFA'
+ include 'COMMON.FFIELD'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SETUP'
C NOTE THAT FILENAMES are FIXED, CURRENTLY!!
open(iodfa, file = 'dist_dfa.dat', status = 'old', err=33)
goto 34
33 write(iout,'(a)') 'Error opening dist_dfa.dat file'
+ print *,'Error opening dist_dfa.dat file'
stop
34 continue
write(iout,'(a)') 'dist_dfa.dat is opened!'
C PHI RESTRAINTS
open(iodfa, file='phi_dfa.dat',status='old',err=35)
goto 36
- 35 write(iout,'(a)') 'Error opening dist_dfa.dat file'
+ 35 write(iout,'(a)') 'Error opening phi_dfa.dat file'
+ print *, 'Error opening phi_dfa.dat file'
stop
36 continue
C THETA RESTRAINTS
open(iodfa, file='theta_dfa.dat',status='old',err=41)
goto 42
- 41 write(iout,'(a)') 'Error opening dist_dfa.dat file'
+ 41 write(iout,'(a)') 'Error opening theta_dfa.dat file'
+ print *,'Error opening theta_dfa.dat file'
stop
42 continue
write(iout,'(a)') 'theta_dfa.dat is opened!'
open(iodfa,file='nei_dfa.dat',status='old',err=37)
goto 38
37 write(iout,'(a)') 'Error opening nei_dfa.dat file'
+ print *,'Error opening nei_dfa.dat file'
stop
38 continue
write(iout,'(a)') 'nei_dfa.dat is opened!'
close(iodfa)
C END OF NEIGHBORING CA
+C init parallel
+C BETA is not parallel !
+#ifdef MPI
+ if (wdfa_beta.ne.0.0 .and. nfgtasks.gt.1) then
+ write (iout,*) "ERRROR dfa_beta works only for FGPROCS=1"
+ print *,"ERRROR dfa_beta works only for FGPROCS=1"
+ stop
+ endif
+ call int_bounds(idfadis,idfadis_start,idfadis_end)
+ call int_bounds(idfaphi,idfaphi_start,idfaphi_end)
+ call int_bounds(idfathe,idfathe_start,idfathe_end)
+ call int_bounds(idfanei,idfanei_start,idfanei_end)
+ if (me.eq.king .or. .not. out1file)
+ & write (iout,*) "DFA MPI ",
+ & "idfadis ",idfadis,idfadis_start,idfadis_end,
+ & "idfaphi ",idfaphi,idfaphi_start,idfaphi_end,
+ & "idfathe ",idfathe,idfathe_start,idfathe_end,
+ & "idfanei ",idfanei,idfanei_start,idfanei_end
+#else
+ idfadis_start=1
+ idfadis_end=idfadis
+ idfaphi_start=1
+ idfaphi_end=idfaphi
+ idfathe_start=1
+ idfathe_end=idfathe
+ idfanei_start=1
+ idfanei_end=idfanei
+#endif
+
+
C READ BETA RESTRAINT
+ if (wdfa_beta.eq.0.0) return
open(iodfa, file='beta_dfa.dat',status='old',err=39)
goto 40
39 write(iout,'(a)') 'Error opening beta_dfa.dat file'
+ print *,'Error opening beta_dfa.dat file'
stop
40 continue
write(iout,'(a)') 'beta_dfa.dat is opened!'
close(iodfa)
C END OF BETA RESTRAINT
+
return
END
edfadis=0
gdfad=0.0d0
- do i=1, idfadis
+ do i=idfadis_start,idfadis_end
iatm1=idislis(1,i)+ishiftca
iatm2=idislis(2,i)+ishiftca
gdfat(:,:) = 0.0d0
C START OF PHI ANGLE
- do i=1, idfaphi
+ do i=idfaphi_start,idfaphi_end
aphi = 0.0d0
do iii=1,5
C END OF PHI RESTRAINT
C START OF THETA ANGLE
- do i=1, idfathe
+ do i=idfathe_start,idfathe_end
athe = 0.0d0
do iii=1,5
c print*, 's1:', s1(:)
c print*, 's2:', s2(:)
- do i=1, idfanei
+ do i=idfanei_start,idfanei_end
kshnum=kshell(i)
n1atom=ineilis(i)+ishiftca
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
projects / unres.git / blobdiff