! Maximum number of SC local term fitting function coefficiants
integer,parameter :: maxsccoef=65
! Maximum number of local shielding effectors
- integer,parameter :: maxcontsshi=50
+! integer,parameter :: maxcontsshi=50
!-----------------------------------------------------------------------------
! commom.calc common/calc/
!-----------------------------------------------------------------------------
#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
! shielding effect varibles for MPI
-! real(kind=8) fac_shieldbuf(maxres),
-! & grad_shield_locbuf(3,maxcontsshi,-1:maxres),
-! & grad_shield_sidebuf(3,maxcontsshi,-1:maxres),
-! & grad_shieldbuf(3,-1:maxres)
-! integer ishield_listbuf(maxres),
-! &shield_listbuf(maxcontsshi,maxres)
+ real(kind=8) fac_shieldbuf(nres), &
+ grad_shield_locbuf(3,maxcontsshi,-1:nres), &
+ grad_shield_sidebuf(3,maxcontsshi,-1:nres), &
+ grad_shieldbuf(3,-1:nres)
+ integer ishield_listbuf(nres), &
+ shield_listbuf(maxcontsshi,nres),k,j,i
! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
! & " nfgtasks",nfgtasks
if (shield_mode.eq.2) then
call set_shield_fac2
endif
+ if (nfgtasks.gt.1) then
+ call MPI_Allgatherv(fac_shield(ivec_start), &
+ ivec_count(fg_rank1), &
+ MPI_DOUBLE_PRECISION,fac_shieldbuf(1),ivec_count(0), &
+ ivec_displ(0), &
+ MPI_DOUBLE_PRECISION,FG_COMM,IERROR)
+ call MPI_Allgatherv(shield_list(1,ivec_start), &
+ ivec_count(fg_rank1), &
+ MPI_I50,shield_listbuf(1,1),ivec_count(0), &
+ ivec_displ(0), &
+ MPI_I50,FG_COMM,IERROR)
+ call MPI_Allgatherv(ishield_list(ivec_start), &
+ ivec_count(fg_rank1), &
+ MPI_INTEGER,ishield_listbuf(1),ivec_count(0), &
+ ivec_displ(0), &
+ MPI_INTEGER,FG_COMM,IERROR)
+ call MPI_Allgatherv(grad_shield(1,ivec_start), &
+ ivec_count(fg_rank1), &
+ MPI_UYZ,grad_shieldbuf(1,1),ivec_count(0), &
+ ivec_displ(0), &
+ MPI_UYZ,FG_COMM,IERROR)
+ call MPI_Allgatherv(grad_shield_side(1,1,ivec_start), &
+ ivec_count(fg_rank1), &
+ MPI_SHI,grad_shield_sidebuf(1,1,1),ivec_count(0), &
+ ivec_displ(0), &
+ MPI_SHI,FG_COMM,IERROR)
+ call MPI_Allgatherv(grad_shield_loc(1,1,ivec_start), &
+ ivec_count(fg_rank1), &
+ MPI_SHI,grad_shield_locbuf(1,1,1),ivec_count(0), &
+ ivec_displ(0), &
+ MPI_SHI,FG_COMM,IERROR)
+ do i=1,nres
+ fac_shield(i)=fac_shieldbuf(i)
+ ishield_list(i)=ishield_listbuf(i)
+ do j=1,3
+ grad_shield(j,i)=grad_shieldbuf(j,i)
+ enddo !j
+ do j=1,ishield_list(i)
+ shield_list(j,i)=shield_listbuf(j,i)
+ do k=1,3
+ grad_shield_loc(k,j,i)=grad_shield_locbuf(k,j,i)
+ grad_shield_side(k,j,i)=grad_shield_sidebuf(k,j,i)
+ enddo !k
+ enddo !j
+ enddo !i
+ endif
+
+
+
+
! print *,"AFTER EGB",ipot,evdw
!mc
!mc Sep-06: egb takes care of dynamic ss bonds too
#ifdef TIMING
time_vec=time_vec+MPI_Wtime()-time01
#endif
+
+
+
+
! print *,"Processor",myrank," left VEC_AND_DERIV"
if (ipot.lt.6) then
#ifdef SPLITELE
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 &
.or. itype(i+3,1).eq.ntyp1 &
.or. itype(i+4,1).eq.ntyp1) cycle
+! print *,"before2",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,i)
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) &
call eturn4(i,eello_turn4)
+! print *,"before",i,i+3, gshieldc_t4(2,i+3),gshieldc_t4(2,i)
num_cont_hb(i)=num_conti
enddo ! i
!
integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield
real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,&
rlocshield
-
+
j=i+3
+! if (j.ne.20) return
+! print *,i,j,gshieldc_t4(2,j),gshieldc_t4(2,j+1)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! Fourth-order contributions
iresshield=shield_list(ilist,i)
do k=1,3
rlocshield=grad_shield_side(k,ilist,i)*eello_t4/fac_shield(i)
+! print *,"rlocshield",rlocshield,grad_shield_side(k,ilist,i),iresshield
gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ &
rlocshield &
+grad_shield_loc(k,ilist,i)*eello_t4/fac_shield(i)
do ilist=1,ishield_list(j)
iresshield=shield_list(ilist,j)
do k=1,3
+! print *,"rlocshieldj",j,rlocshield,grad_shield_side(k,ilist,j),iresshield
rlocshield=grad_shield_side(k,ilist,j)*eello_t4/fac_shield(j)
gshieldx_t4(k,iresshield)=gshieldx_t4(k,iresshield)+ &
rlocshield &
+grad_shield_loc(k,ilist,j)*eello_t4/fac_shield(j)
gshieldc_t4(k,iresshield-1)=gshieldc_t4(k,iresshield-1) &
+rlocshield
+! print *,"after", gshieldc_t4(k,iresshield-1),iresshield-1,gshieldc_t4(k,iresshield)
enddo
enddo
-
do k=1,3
gshieldc_t4(k,i)=gshieldc_t4(k,i)+ &
grad_shield(k,i)*eello_t4/fac_shield(i)
grad_shield(k,i)*eello_t4/fac_shield(i)
gshieldc_t4(k,j-1)=gshieldc_t4(k,j-1)+ &
grad_shield(k,j)*eello_t4/fac_shield(j)
+! print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1)
enddo
endif
call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
+! if (j.lt.nres-1) then
gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) &
*fac_shield(i)*fac_shield(j) &
*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
-
+! endif
a_temp(1,1)=aggj1(l,1)
a_temp(1,2)=aggj1(l,2)
call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
s3=0.5d0*(pizda(1,1)+pizda(2,2))
! write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
+! if (j.lt.nres-1) then
+! print *,"juest before",j1, gcorr4_turn(l,j1)
gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) &
*fac_shield(i)*fac_shield(j) &
*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
-
+! if (shield_mode.gt.0) then
+! print *,"juest after",j1, gcorr4_turn(l,j1),gshieldc_t4(k,j1),gshieldc_loc_t4(k,j1),gel_loc_turn4(i+2)
+! else
+! print *,"juest after",j1, gcorr4_turn(l,j1),gel_loc_turn4(i+2)
+! endif
+! endif
enddo
gshieldc_t4(3,i)=gshieldc_t4(3,i)+ &
ssgradlipi*eello_t4/4.0d0*lipscale
wscpho*gvdwc_scpho(j,i)+ &
wpeppho*gvdwc_peppho(j,i)
-
+
+0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)&
+wvdwpsb*gvdwpsb1(j,i))&
+wbond_nucl*gradb_nucl(j,i)+wsbloc*gsbloc(j,i)
-
+! if (i.eq.21) then
+! print *,"in sum",gradc(j,i,icg),wturn4*gcorr4_turn(j,i),&
+! wturn4*gshieldc_t4(j,i), &
+! wturn4*gshieldc_loc_t4(j,i)
+! endif
! if ((i.le.2).and.(i.ge.1))
! print *,gradc(j,i,icg),&
! gradbufc(j,i),welec*gelc(j,i), &
! if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
enddo
- enddo
+ enddo
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "gloc before adding corr"
do i=1,4*nres
write (iout,*) i,gloc(i,icg)
enddo
#endif
+!#undef DEBUG
#ifdef MPI
if (nfgtasks.gt.1) then
do j=1,3
endif
endif
endif
-!el#define DEBUG
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "gradc gradx gloc"
do i=1,nres
i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg)
enddo
#endif
-!el#undef DEBUG
+!#undef DEBUG
#ifdef TIMING
time_sumgradient=time_sumgradient+MPI_Wtime()-time01
#endif
c(j,k)=c(j,k)+aincr
c(j,k+nres)=c(j,k+nres)+aincr
enddo
+ call zerograd
call etotal(energia1)
etot1=energia1(0)
ggg(j)=(etot1-etot)/aincr
do j=1,3
c(j,i+nres)=c(j,i+nres)+aincr
dc(j,i+nres)=dc(j,i+nres)+aincr
+ call zerograd
call etotal(energia1)
etot1=energia1(0)
ggg(j+3)=(etot1-etot)/aincr
write (iout,*) "split_ene ",split_ene
call flush(iout)
if (.not.split_ene) then
+ call zerograd
call etotal(energia)
etot=energia(0)
call cartgrad
dc(j,i+nres)=c(j,i+nres)-c(j,i)
call int_from_cart1(.false.)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot1=energia1(0)
write (iout,*) "ij",i,j," etot1",etot1
dc(j,i+nres)=c(j,i+nres)-c(j,i)
call int_from_cart1(.false.)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot2=energia1(0)
write (iout,*) "ij",i,j," etot2",etot2
dc(j,i+nres)=c(j,i+nres)-c(j,i)
call int_from_cart1(.false.)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot1=energia1(0)
else
dc(j,i+nres)=c(j,i+nres)-c(j,i)
call int_from_cart1(.false.)
if (.not.split_ene) then
- call etotal(energia1)
+ call zerograd
+ call etotal(energia1)
etot2=energia1(0)
ggg(j+3)=(etot1-etot2)/(2*aincr)
else
do i=1,nres
do j=1,3
grad_s(j,i)=gcart(j,i)
+! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i)
+
! if (i.le.2) print *,"tu?!",gcart(j,i),grad_s(j,i),gxcart(j,i)
grad_s(j+3,i)=gxcart(j,i)
enddo
do i=1,nres
do j=1,3
grad_s(j,i)=gcart(j,i)
+! if (i.eq.21) print *,"PRZEKAZANIE",gcart(j,i)
grad_s(j+3,i)=gxcart(j,i)
enddo
enddo
#endif
! call int_from_cart1(.false.)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot1=energia1(0)
! call enerprint(energia1)
call chainbuild_cart
! call int_from_cart1(.false.)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot2=energia1(0)
ggg(j)=(etot1-etot2)/(2*aincr)
! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
! write (iout,*)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot1=energia1(0)
else
! write (iout,*) "dxnormnormsafe",dsqrt(
! & dxnorm_safe(1)**2+dxnorm_safe(2)**2+dxnorm_safe(3)**2)
if (.not.split_ene) then
+ call zerograd
call etotal(energia1)
etot2=energia1(0)
ggg(j+3)=(etot1-etot2)/(2*aincr)
call sum_gradient
#ifdef TIMING
#endif
+!#define DEBUG
!el write (iout,*) "After sum_gradient"
#ifdef DEBUG
!el write (iout,*) "After sum_gradient"
write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
enddo
#endif
+!#undef DEBUG
! If performing constraint dynamics, add the gradients of the constraint energy
if(usampl.and.totT.gt.eq_time) then
do i=1,nct
#endif
! call checkintcartgrad
! write(iout,*) 'calling int_to_cart'
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "gcart, gxcart, gloc before int_to_cart"
#endif
(gxcart(j,i),j=1,3)
enddo
#endif
+!#undef DEBUG
#ifdef CARGRAD
#ifdef DEBUG
write (iout,*) "CARGRAD"
dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+ &
cost1*dc_norm(j,i-2))/ &
vbld(i-1)
- domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)
+ domicron(j,1,1,i)=-1.0/sint1*dcosomicron(j,1,1,i)
dcosomicron(j,1,2,i)=-(dc_norm(j,i-2) &
+cost1*(dc_norm(j,i-1+nres)))/ &
vbld(i-1+nres)
- domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)
+ domicron(j,1,2,i)=-1.0/sint1*dcosomicron(j,1,2,i)
!C Calculate derivative over second omicron Sci-1,Cai-1 Cai
!C Looks messy but better than if in loop
dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres) &
+cost2*dc_norm(j,i-1))/ &
vbld(i)
- domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
+ domicron(j,2,1,i)=-1.0/sint2*dcosomicron(j,2,1,i)
dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) &
+cost2*(-dc_norm(j,i-1+nres)))/ &
vbld(i-1+nres)
! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres)
- domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)
+ domicron(j,2,2,i)=-1.0/sint2*dcosomicron(j,2,2,i)
enddo
endif
enddo
dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
dc_norm(j,i-3))/vbld(i-2)
- dphi(j,1,i)=-1/sing*dcosphi(j,1,i)
+ dphi(j,1,i)=-1.0/sing*dcosphi(j,1,i)
dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2* &
dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4* &
dcostheta(j,1,i)
- dphi(j,2,i)=-1/sing*dcosphi(j,2,i)
+ dphi(j,2,i)=-1.0/sing*dcosphi(j,2,i)
dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4* &
dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp* &
dc_norm(j,i-1))/vbld(i)
- dphi(j,3,i)=-1/sing*dcosphi(j,3,i)
+ dphi(j,3,i)=-1.0/sing*dcosphi(j,3,i)
+!#define DEBUG
+#ifdef DEBUG
+ write(iout,*) "just after",dphi(j,3,i),sing,dcosphi(j,3,i)
+#endif
+!#undef DEBUG
endif
enddo
endif
call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),&
MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,&
king,FG_COMM,IERROR)
+!#define DEBUG
#ifdef DEBUG
!d write (iout,*) "Gather dphi"
!d call flush(iout)
write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
enddo
#endif
+!#undef DEBUG
call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),&
MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,&
king,FG_COMM,IERROR)
#endif
endif
#endif
+!#define DEBUG
#ifdef DEBUG
write (iout,*) "dtheta after gather"
do i=1,nres
write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
enddo
#endif
+!#undef DEBUG
return
end subroutine intcartderiv
!-----------------------------------------------------------------------------
enddo
!C now sscale fraction
sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
-!C print *,buff_shield,"buff"
+! print *,buff_shield,"buff",sh_frac_dist
!C now sscale
if (sh_frac_dist.le.0.0) cycle
!C print *,ishield_list(i),i
long=long_r_sidechain(itype(k,1))
costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
sinthet=short/dist_pep_side*costhet
+! print *,"SORT",short,long,sinthet,costhet
!C now costhet_grad
!C costhet=0.6d0
!C sinthet=0.8
enddo
!C print *,sinphi,sinthet
VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) &
- & /VSolvSphere_div
+ /VSolvSphere_div
!C & *wshield
!C now the gradient...
do j=1,3
sinphi/sinthet*costhet*costhet_grad(j)&
+sinthet/sinphi*cosphi*cosphi_grad_long(j))) &
)*wshield
-
+! print *, 1.0d0/(-dsqrt(1.0d0-sinphi*sinthet)),&
+! sinphi/sinthet,&
+! +sinthet/sinphi,"HERE"
grad_shield_loc(j,ishield_list(i),i)= &
scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*&
(1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(&
sinthet/sinphi*cosphi*cosphi_grad_loc(j)&
))&
*wshield
+! print *,grad_shield_loc(j,ishield_list(i),i)
enddo
VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
enddo
fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
-!C write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i)
+! write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i)
enddo
return
end subroutine set_shield_fac2
projects / unres4.git / blobdiff