grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres)
!-----------------------------NUCLEIC GRADIENT
real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl, &
- gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc
+ gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc,gsbloc,&
+ gsblocx,gradcorr_nucl,gradxorr_nucl,gradcorr3_nucl,gradxorr3_nucl
! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,&
gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres)
real(kind=8),dimension(:,:,:,:),allocatable :: uygrad,uzgrad !(3,3,2,maxres)
!-----------------------------------------------------------------------------
! common /przechowalnia/
- real(kind=8),dimension(:,:,:),allocatable :: zapas !(max_dim,maxconts,max_fg_procs)
+ real(kind=8),dimension(:,:,:),allocatable :: zapas
+ real(kind=8),dimension(:,:,:,:),allocatable ::zapas2 !(max_dim,maxconts,max_fg_procs)
real(kind=8),dimension(:,:,:),allocatable :: fromto !(3,3,maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
etube=0.0d0
endif
!--------------------------------------------------------
+ print *,"before",ees,evdw1,ecorr
call ebond_nucl(estr_nucl)
call ebend_nucl(ebe_nucl)
call etor_nucl(etors_nucl)
call esb_gb(evdwsb,eelsb)
-! call multibody_hb(ecorr,ecorr3,n_corr,n_corr1)
call epp_nucl_sub(evdwpp,eespp)
call epsb(evdwpsb,eelpsb)
+ call esb(esbloc)
+ call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
print *,"after ebend", ebe_nucl
#ifdef TIMING
!
! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
!
+ print *,"iatel_s,iatel_e,",iatel_s,iatel_e
do i=iatel_s,iatel_e
if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle
dxi=dc(1,i)
iteli=itel(i)
itelj=itel(j)
if (j.eq.i+2 .and. itelj.eq.2) iteli=2
- aaa=app_nucl(iteli,itelj)
- bbb=bpp_nucl(iteli,itelj)
- ael6i=ael6_nucl(iteli,itelj)
- ael3i=ael3_nucl(iteli,itelj)
+ aaa=app(iteli,itelj)
+ bbb=bpp(iteli,itelj)
+ ael6i=ael6(iteli,itelj)
+ ael3i=ael3(iteli,itelj)
dxj=dc(1,j)
dyj=dc(2,j)
dzj=dc(3,j)
allocate(gelsbx(3,-1:nres))
allocate(gvdwsbx(3,-1:nres))
allocate(gvdwsbc(3,-1:nres))
+ allocate(gsbloc(3,-1:nres))
+ allocate(gsblocx(3,-1:nres))
+ allocate(gradcorr_nucl(3,-1:nres))
+ allocate(gradxorr_nucl(3,-1:nres))
+ allocate(gradcorr3_nucl(3,-1:nres))
+ allocate(gradxorr3_nucl(3,-1:nres))
+
!(3,maxres)
allocate(grad_shield_side(3,50,nres))
allocate(grad_shield_loc(3,50,nres))
subroutine esb_gb(evdwsb,eelsb)
use comm_locel
use calc_data_nucl
- integer :: iint,itypi,itypi1,itypj,subchap
+ integer :: iint,itypi,itypi1,itypj,subchap,num_conti2
real(kind=8) :: xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0iji,evdwsb,eelsb,ecorr,eelij
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
! print *,"iastsc_nucl",iatsc_s_nucl,iatsc_e_nucl
do i=iatsc_s_nucl,iatsc_e_nucl
num_conti=0
+ num_conti2=0
itypi=itype(i,2)
! PRINT *,"I=",i,itypi
if (itypi.eq.ntyp1_molec(2)) cycle
gg(3)=zj*fac
!C Calculate angular part of the gradient.
call sc_grad_nucl
- call eelsbij(eelij)
+ call eelsbij(eelij,num_conti2)
if (energy_dec .and. &
(j.eq.i+1.or.j.eq.nres-i+1.or.j.eq.nres-i.or.j.eq.nres-i+2)) &
write (istat,'(e14.5)') evdwij
eelsb=eelsb+eelij
enddo ! j
enddo ! iint
- num_cont_hb(i)=num_conti
+ num_cont_hb(i)=num_conti2
enddo ! i
!c write (iout,*) "Number of loop steps in EGB:",ind
!cccc energy_dec=.false.
return
end subroutine esb_gb
!-------------------------------------------------------------------------------
- subroutine eelsbij(eesij)
+ subroutine eelsbij(eesij,num_conti2)
use comm_locel
use calc_data_nucl
real(kind=8),dimension(3) :: ggg,gggp,gggm,dcosb,dcosg
real(kind=8),dimension(3,3) :: erder,uryg,urzg,vryg,vrzg
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dist_temp, dist_init,rlocshield,fracinbuf
- integer xshift,yshift,zshift,ilist,iresshield
+ integer xshift,yshift,zshift,ilist,iresshield,num_conti2
!c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
real(kind=8) scal_el /0.5d0/
gelsbc(k,j)=gelsbc(k,j)+ggg(k)
gelsbc(k,i)=gelsbc(k,i)-ggg(k)
enddo
- IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and. j.gt.i+1 .and.&
+! IF ( (wcorr_nucl.gt.0.0d0.or.wcorr3_nucl.gt.0.0d0) .and.
+ IF ( j.gt.i+1 .and.&
num_conti.le.maxconts) THEN
!C
!C Calculate the contact function. The ith column of the array JCONT will
!c write (2,*) "fcont",fcont
if (fcont.gt.0.0D0) then
num_conti=num_conti+1
+ num_conti2=num_conti2+1
+
if (num_conti.gt.maxconts) then
write (iout,*) 'WARNING - max. # of contacts exceeded;',&
' will skip next contacts for this conf.'
enddo
return
end subroutine sc_grad_nucl
+!-----------------------------------------------------------------------
+ subroutine esb(esbloc)
+!C Calculate the local energy of a side chain and its derivatives in the
+!C corresponding virtual-bond valence angles THETA and the spherical angles
+!C ALPHA and OMEGA derived from AM1 all-atom calculations.
+!C added by Urszula Kozlowska. 07/11/2007
+!C
+ real(kind=8),dimension(3):: x_prime,y_prime,z_prime
+ real(kind=8),dimension(9):: x
+ real(kind=8) :: sumene,dsc_i,dp2_i,xx,yy,zz,sumene1, &
+ sumene2,sumene3,sumene4,s1,s1_6,s2,s2_6,&
+ de_dxx,de_dyy,de_dzz,de_dt,s1_t,s1_6_t,s2_t,s2_6_t
+ real(kind=8),dimension(3):: dXX_Ci1,dYY_Ci1,dZZ_Ci1,dXX_Ci,&
+ dYY_Ci,dZZ_Ci,dXX_XYZ,dYY_XYZ,dZZ_XYZ,dt_dCi,dt_dCi1
+ real(kind=8) :: esbloc,delta,cosfac2,cosfac,sinfac2,sinfac,de_dtt,&
+ cossc,cossc1,cosfac2xx,sinfac2yy,pom1,pom
+ integer::it,nlobit,i,j,k
+! common /sccalc/ time11,time12,time112,theti,it,nlobit
+ delta=0.02d0*pi
+ esbloc=0.0D0
+ do i=loc_start_nucl,loc_end_nucl
+ if (itype(i,2).eq.ntyp1_molec(2)) cycle
+ costtab(i+1) =dcos(theta(i+1))
+ sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+ cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+ sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+ cosfac2=0.5d0/(1.0d0+costtab(i+1))
+ cosfac=dsqrt(cosfac2)
+ sinfac2=0.5d0/(1.0d0-costtab(i+1))
+ sinfac=dsqrt(sinfac2)
+ it=itype(i,2)
+ if (it.eq.10) goto 1
+
+!c
+!C Compute the axes of tghe local cartesian coordinates system; store in
+!c x_prime, y_prime and z_prime
+!c
+ do j=1,3
+ x_prime(j) = 0.00
+ y_prime(j) = 0.00
+ z_prime(j) = 0.00
+ enddo
+!C write(2,*) "dc_norm", dc_norm(1,i+nres),dc_norm(2,i+nres),
+!C & dc_norm(3,i+nres)
+ do j = 1,3
+ x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+ enddo
+ do j = 1,3
+ z_prime(j) = -uz(j,i-1)
+ enddo
+
+ xx=0.0d0
+ yy=0.0d0
+ zz=0.0d0
+ do j = 1,3
+ xx = xx + x_prime(j)*dc_norm(j,i+nres)
+ yy = yy + y_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ enddo
+
+ xxtab(i)=xx
+ yytab(i)=yy
+ zztab(i)=zz
+ it=itype(i,2)
+ do j = 1,9
+ x(j) = sc_parmin_nucl(j,it)
+ enddo
+#ifdef CHECK_COORD
+!Cc diagnostics - remove later
+ xx1 = dcos(alph(2))
+ yy1 = dsin(alph(2))*dcos(omeg(2))
+ zz1 = -dsin(alph(2))*dsin(omeg(2))
+ write(2,'(3f8.1,3f9.3,1x,3f9.3)') &
+ alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,&
+ xx1,yy1,zz1
+!C," --- ", xx_w,yy_w,zz_w
+!c end diagnostics
+#endif
+ sumene = enesc_nucl(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ esbloc = esbloc + sumene
+ if (energy_dec) write(iout,*) "i",i," esbloc",sumene,esbloc,xx,yy,zz
+ if (energy_dec) write(iout,*) "x",(x(k),k=1,9)
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+!C
+!C This section to check the numerical derivatives of the energy of ith side
+!C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert
+!C #define DEBUG in the code to turn it on.
+!C
+ write (2,*) "sumene =",sumene
+ aincr=1.0d-7
+ xxsave=xx
+ xx=xx+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dxx_num=(sumenep-sumene)/aincr
+ xx=xxsave
+ write (2,*) "xx+ sumene from enesc=",sumenep,sumene
+ yysave=yy
+ yy=yy+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dyy_num=(sumenep-sumene)/aincr
+ yy=yysave
+ write (2,*) "yy+ sumene from enesc=",sumenep,sumene
+ zzsave=zz
+ zz=zz+aincr
+ write (2,*) xx,yy,zz
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dzz_num=(sumenep-sumene)/aincr
+ zz=zzsave
+ write (2,*) "zz+ sumene from enesc=",sumenep,sumene
+ costsave=cost2tab(i+1)
+ sintsave=sint2tab(i+1)
+ cost2tab(i+1)=dcos(0.5d0*(theta(i+1)+aincr))
+ sint2tab(i+1)=dsin(0.5d0*(theta(i+1)+aincr))
+ sumenep = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1))
+ de_dt_num=(sumenep-sumene)/aincr
+ write (2,*) " t+ sumene from enesc=",sumenep,sumene
+ cost2tab(i+1)=costsave
+ sint2tab(i+1)=sintsave
+!C End of diagnostics section.
+#endif
+!C
+!C Compute the gradient of esc
+!C
+ de_dxx=x(1)+2*x(4)*xx+x(7)*zz+x(8)*yy
+ de_dyy=x(2)+2*x(5)*yy+x(8)*xx+x(9)*zz
+ de_dzz=x(3)+2*x(6)*zz+x(7)*xx+x(9)*yy
+ de_dtt=0.0d0
+#ifdef DEBUG
+ write (2,*) "x",(x(k),k=1,9)
+ write (2,*) "xx",xx," yy",yy," zz",zz
+ write (2,*) "de_xx ",de_xx," de_yy ",de_yy,&
+ " de_zz ",de_zz," de_tt ",de_tt
+ write (2,*) "de_xx_num",de_dxx_num," de_yy_num",de_dyy_num,&
+ " de_zz_num",de_dzz_num," de_dt_num",de_dt_num
+#endif
+!C
+ cossc=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ cossc1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ cosfac2xx=cosfac2*xx
+ sinfac2yy=sinfac2*yy
+ do k = 1,3
+ dt_dCi(k) = -(dc_norm(k,i-1)+costtab(i+1)*dc_norm(k,i))*&
+ vbld_inv(i+1)
+ dt_dCi1(k)= -(dc_norm(k,i)+costtab(i+1)*dc_norm(k,i-1))*&
+ vbld_inv(i)
+ pom=(dC_norm(k,i+nres)-cossc*dC_norm(k,i))*vbld_inv(i+1)
+ pom1=(dC_norm(k,i+nres)-cossc1*dC_norm(k,i-1))*vbld_inv(i)
+!c write (iout,*) "i",i," k",k," pom",pom," pom1",pom1,
+!c & " dt_dCi",dt_dCi(k)," dt_dCi1",dt_dCi1(k)
+!c write (iout,*) "dC_norm",(dC_norm(j,i),j=1,3),
+!c & (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
+ dXX_Ci(k)=pom*cosfac-dt_dCi(k)*cosfac2xx
+ dXX_Ci1(k)=-pom1*cosfac-dt_dCi1(k)*cosfac2xx
+ dYY_Ci(k)=pom*sinfac+dt_dCi(k)*sinfac2yy
+ dYY_Ci1(k)=pom1*sinfac+dt_dCi1(k)*sinfac2yy
+ dZZ_Ci1(k)=0.0d0
+ dZZ_Ci(k)=0.0d0
+ do j=1,3
+ dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
+ dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+ enddo
+
+ dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
+ dYY_XYZ(k)=vbld_inv(i+nres)*(y_prime(k)-yy*dC_norm(k,i+nres))
+ dZZ_XYZ(k)=vbld_inv(i+nres)*(z_prime(k)-zz*dC_norm(k,i+nres))
+!c
+ dt_dCi(k) = -dt_dCi(k)/sinttab(i+1)
+ dt_dCi1(k)= -dt_dCi1(k)/sinttab(i+1)
+ enddo
+
+ do k=1,3
+ dXX_Ctab(k,i)=dXX_Ci(k)
+ dXX_C1tab(k,i)=dXX_Ci1(k)
+ dYY_Ctab(k,i)=dYY_Ci(k)
+ dYY_C1tab(k,i)=dYY_Ci1(k)
+ dZZ_Ctab(k,i)=dZZ_Ci(k)
+ dZZ_C1tab(k,i)=dZZ_Ci1(k)
+ dXX_XYZtab(k,i)=dXX_XYZ(k)
+ dYY_XYZtab(k,i)=dYY_XYZ(k)
+ dZZ_XYZtab(k,i)=dZZ_XYZ(k)
+ enddo
+ do k = 1,3
+!c write (iout,*) "k",k," dxx_ci1",dxx_ci1(k)," dyy_ci1",
+!c & dyy_ci1(k)," dzz_ci1",dzz_ci1(k)
+!c write (iout,*) "k",k," dxx_ci",dxx_ci(k)," dyy_ci",
+!c & dyy_ci(k)," dzz_ci",dzz_ci(k)
+!c write (iout,*) "k",k," dt_dci",dt_dci(k)," dt_dci",
+!c & dt_dci(k)
+!c write (iout,*) "k",k," dxx_XYZ",dxx_XYZ(k)," dyy_XYZ",
+!c & dyy_XYZ(k)," dzz_XYZ",dzz_XYZ(k)
+ gsbloc(k,i-1)=gsbloc(k,i-1)+de_dxx*dxx_ci1(k) &
+ +de_dyy*dyy_ci1(k)+de_dzz*dzz_ci1(k)+de_dt*dt_dCi1(k)
+ gsbloc(k,i)=gsbloc(k,i)+de_dxx*dxx_Ci(k) &
+ +de_dyy*dyy_Ci(k)+de_dzz*dzz_Ci(k)+de_dt*dt_dCi(k)
+ gsblocx(k,i)= de_dxx*dxx_XYZ(k)&
+ +de_dyy*dyy_XYZ(k)+de_dzz*dzz_XYZ(k)
+ enddo
+!c write(iout,*) "ENERGY GRAD = ", (gsbloc(k,i-1),k=1,3),
+!c & (gsbloc(k,i),k=1,3),(gsblocx(k,i),k=1,3)
+
+!C to check gradient call subroutine check_grad
+
+ 1 continue
+ enddo
+ return
+ end subroutine esb
+!=-------------------------------------------------------
+ real(kind=8) function enesc_nucl(x,xx,yy,zz,cost2,sint2)
+! implicit none
+ real(kind=8),dimension(9):: x(9)
+ real(kind=8) :: xx,yy,zz,cost2,sint2,sumene1,sumene2, &
+ sumene3,sumene4,sumene,dsc_i,dp2_i,dscp1,dscp2,s1,s1_6,s2,s2_6
+ integer i
+!c write (2,*) "enesc"
+!c write (2,*) "x",(x(i),i=1,9)
+!c write(2,*)"xx",xx," yy",yy," zz",zz," cost2",cost2," sint2",sint2
+ sumene=x(1)*xx+x(2)*yy+x(3)*zz+x(4)*xx**2 &
+ + x(5)*yy**2+x(6)*zz**2+x(7)*xx*zz+x(8)*xx*yy &
+ + x(9)*yy*zz
+ enesc_nucl=sumene
+ return
+ end function enesc_nucl
+!-----------------------------------------------------------------------------
+ subroutine multibody_hb_nucl(ecorr,ecorr3,n_corr,n_corr1)
+#ifdef MPI
+ include 'mpif.h'
+ integer,parameter :: max_cont=2000
+ integer,parameter:: max_dim=2*(8*3+6)
+ integer, parameter :: msglen1=max_cont*max_dim
+ integer,parameter :: msglen2=2*msglen1
+ integer source,CorrelType,CorrelID,Error
+ real(kind=8) :: buffer(max_cont,max_dim)
+ integer status(MPI_STATUS_SIZE)
+ integer :: ierror,nbytes
+#endif
+ real(kind=8),dimension(3):: gx(3),gx1(3)
+ real(kind=8) :: time00
+ logical lprn,ldone
+ integer i,j,i1,j1,jj,kk,num_conti,num_conti1,nn
+ real(kind=8) ecorr,ecorr3
+ integer :: n_corr,n_corr1,mm,msglen
+!C Set lprn=.true. for debugging
+ lprn=.false.
+ n_corr=0
+ n_corr1=0
+#ifdef MPI
+ if(.not.allocated(zapas2)) allocate(zapas2(3,maxconts,nres,8))
+
+ if (nfgtasks.le.1) goto 30
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt,nct-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+!C Caution! Following code assumes that electrostatic interactions concerning
+!C a given atom are split among at most two processors!
+ CorrelType=477
+ CorrelID=fg_rank+1
+ ldone=.false.
+ do i=1,max_cont
+ do j=1,max_dim
+ buffer(i,j)=0.0D0
+ enddo
+ enddo
+ mm=mod(fg_rank,2)
+!c write (*,*) 'MyRank',MyRank,' mm',mm
+ if (mm) 20,20,10
+ 10 continue
+!c write (*,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.gt.0) then
+!C Send correlation contributions to the preceding processor
+ msglen=msglen1
+ nn=num_cont_hb(iatel_s_nucl)
+ call pack_buffer(max_cont,max_dim,iatel_s,0,buffer)
+!c write (*,*) 'The BUFFER array:'
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,30)
+!c enddo
+ if (ielstart_nucl(iatel_s_nucl).gt.iatel_s_nucl+ispp) then
+ msglen=msglen2
+ call pack_buffer(max_cont,max_dim,iatel_s+1,30,buffer)
+!C Clear the contacts of the atom passed to the neighboring processor
+ nn=num_cont_hb(iatel_s_nucl+1)
+!c do i=1,nn
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j+30),j=1,30)
+!c enddo
+ num_cont_hb(iatel_s_nucl)=0
+ endif
+!cd write (iout,*) 'Processor ',fg_rank,MyRank,
+!cd & ' is sending correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen
+!c write (*,*) 'Processor ',fg_rank,MyRank,
+!c & ' is sending correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ call MPI_Send(buffer,msglen,MPI_DOUBLE_PRECISION,fg_rank-1, &
+ CorrelType,FG_COMM,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!cd write (iout,*) 'Processor ',fg_rank,
+!cd & ' has sent correlation contribution to processor',fg_rank-1,
+!cd & ' msglen=',msglen,' CorrelID=',CorrelID
+!c write (*,*) 'Processor ',fg_rank,
+!c & ' has sent correlation contribution to processor',fg_rank-1,
+!c & ' msglen=',msglen,' CorrelID=',CorrelID
+!c msglen=msglen1
+ endif ! (fg_rank.gt.0)
+ if (ldone) goto 30
+ ldone=.true.
+ 20 continue
+!c write (*,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone
+ if (fg_rank.lt.nfgtasks-1) then
+!C Receive correlation contributions from the next processor
+ msglen=msglen1
+ if (ielend_nucl(iatel_e_nucl).lt.nct_molec(2)-1) msglen=msglen2
+!cd write (iout,*) 'Processor',fg_rank,
+!cd & ' is receiving correlation contribution from processor',fg_rank+1,
+!cd & ' msglen=',msglen,' CorrelType=',CorrelType
+!c write (*,*) 'Processor',fg_rank,
+!c &' is receiving correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' CorrelType=',CorrelType
+ time00=MPI_Wtime()
+ nbytes=-1
+ do while (nbytes.le.0)
+ call MPI_Probe(fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ call MPI_Get_count(status,MPI_DOUBLE_PRECISION,nbytes,IERROR)
+ enddo
+!c print *,'Processor',myrank,' msglen',msglen,' nbytes',nbytes
+ call MPI_Recv(buffer,nbytes,MPI_DOUBLE_PRECISION, &
+ fg_rank+1,CorrelType,FG_COMM,status,IERROR)
+ time_sendrecv=time_sendrecv+MPI_Wtime()-time00
+!c write (*,*) 'Processor',fg_rank,
+!c &' has received correlation contribution from processor',fg_rank+1,
+!c & ' msglen=',msglen,' nbytes=',nbytes
+!c write (*,*) 'The received BUFFER array:'
+!c do i=1,max_cont
+!c write (*,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,60)
+!c enddo
+ if (msglen.eq.msglen1) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,0,buffer)
+ else if (msglen.eq.msglen2) then
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl,0,buffer)
+ call unpack_buffer(max_cont,max_dim,iatel_e_nucl+1,30,buffer)
+ else
+ write (iout,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ write (*,*) &
+ 'ERROR!!!! message length changed while processing correlations.'
+ call MPI_Abort(MPI_COMM_WORLD,Error,IERROR)
+ endif ! msglen.eq.msglen1
+ endif ! fg_rank.lt.nfgtasks-1
+ if (ldone) goto 30
+ ldone=.true.
+ goto 10
+ 30 continue
+#endif
+ if (lprn) then
+ write (iout,'(a)') 'Contact function values:'
+ do i=nnt_molec(2),nct_molec(2)-1
+ write (iout,'(2i3,50(1x,i2,f5.2))') &
+ i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), &
+ j=1,num_cont_hb(i))
+ enddo
+ endif
+ ecorr=0.0D0
+ ecorr3=0.0d0
+!C Remove the loop below after debugging !!!
+ do i=nnt_molec(2),nct_molec(2)
+ do j=1,3
+ gradcorr_nucl(j,i)=0.0D0
+ gradxorr_nucl(j,i)=0.0D0
+ gradcorr3_nucl(j,i)=0.0D0
+ gradxorr3_nucl(j,i)=0.0D0
+ enddo
+ enddo
+ print *,"iatsc_s_nucl,iatsc_e_nucl",iatsc_s_nucl,iatsc_e_nucl
+!C Calculate the local-electrostatic correlation terms
+ do i=iatsc_s_nucl,iatsc_e_nucl
+ i1=i+1
+ num_conti=num_cont_hb(i)
+ num_conti1=num_cont_hb(i+1)
+ print *,i,num_conti,num_conti1
+ do jj=1,num_conti
+ j=jcont_hb(jj,i)
+ do kk=1,num_conti1
+ j1=jcont_hb(kk,i1)
+!c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1 .or. j1.eq.j-1) then
+!C
+!C Contacts I-J and (I+1)-(J+1) or (I+1)-(J-1) occur simultaneously.
+!C The system gains extra energy.
+!C Tentative expression & coefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 34 and 35 from Liwo et al., 1998.
+!C
+ ecorr=ecorr+ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') &
+ 'ecorrh',i,j,ehbcorr_nucl(i,j,i+1,j1,jj,kk,0.528D0,0.132D0)
+ n_corr=n_corr+1
+ else if (j1.eq.j) then
+!C
+!C Contacts I-J and I-(J+1) occur simultaneously.
+!C The system loses extra energy.
+!C Tentative expression & c?oefficients; assumed d(stacking)=4.5 A,
+!C parallel dipoles of stacknig bases and sin(mui)sin(muj)/eps/d^3=0.7
+!C Need to implement full formulas 32 from Liwo et al., 1998.
+!C
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i+1,j,jj,kk,0.310D0,-0.155D0)
+ endif
+ enddo ! kk
+ do kk=1,num_conti
+ j1=jcont_hb(kk,i)
+!c write (iout,*) 'ecorr3: i=',i,' j=',j,' i1=',i1,' j1=',j1,
+!c & ' jj=',jj,' kk=',kk
+ if (j1.eq.j+1) then
+!C Contacts I-J and (I+1)-J occur simultaneously.
+!C The system loses extra energy.
+ ecorr3=ecorr3+ehbcorr3_nucl(i,j,i,j+1,jj,kk,0.310D0,-0.155D0)
+ endif ! j1==j+1
+ enddo ! kk
+ enddo ! jj
+ enddo ! i
+ return
+ end subroutine multibody_hb_nucl
+!-----------------------------------------------------------
+ real(kind=8) function ehbcorr_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
+ real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
+ ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+! print *,"ehbcorr_nucl",ekont,ees
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr_nucl=ecorr_nucl+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+&
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr_nucl(ll,j)=gradcorr_nucl(ll,j)+gradlongij
+ gradcorr_nucl(ll,i)=gradcorr_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr_nucl(ll,l)=gradcorr_nucl(ll,l)+gradlongkl
+ gradcorr_nucl(ll,k)=gradcorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,i)=gradxorr_nucl(ll,i)-gradlongij
+ gradxorr_nucl(ll,j)=gradxorr_nucl(ll,j)+gradlongij
+ gradxorr_nucl(ll,k)=gradxorr_nucl(ll,k)-gradlongkl
+ gradxorr_nucl(ll,l)=gradxorr_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr_nucl=ekont*ees
+ return
+ end function ehbcorr_nucl
+!-------------------------------------------------------------------------
+
+ real(kind=8) function ehbcorr3_nucl(i,j,k,l,jj,kk,coeffp,coeffm)
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.DERIV'
+! include 'COMMON.INTERACT'
+! include 'COMMON.CONTACTS'
+ real(kind=8),dimension(3) :: gx,gx1
+ logical :: lprn
+!el local variables
+ integer :: i,j,k,l,jj,kk,ll,ilist,m, iresshield
+ real(kind=8) :: coeffp,coeffm,eij,ekl,ees0pij,ees0pkl,ees0mij,&
+ ees0mkl,ees,coeffpees0pij,coeffmees0mij,&
+ coeffpees0pkl,coeffmees0mkl,gradlongij,gradlongkl, &
+ rlocshield
+
+ lprn=.false.
+ eij=facont_hb(jj,i)
+ ekl=facont_hb(kk,k)
+ ees0pij=ees0p(jj,i)
+ ees0pkl=ees0p(kk,k)
+ ees0mij=ees0m(jj,i)
+ ees0mkl=ees0m(kk,k)
+ ekont=eij*ekl
+ ees=-(coeffp*ees0pij*ees0pkl+coeffm*ees0mij*ees0mkl)
+!cd ees=-(coeffp*ees0pkl+coeffm*ees0mkl)
+!C Following 4 lines for diagnostics.
+!cd ees0pkl=0.0D0
+!cd ees0pij=1.0D0
+!cd ees0mkl=0.0D0
+!cd ees0mij=1.0D0
+!cd write (iout,*)'Contacts have occurred for nucleic bases',
+!cd & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l
+!cd & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees
+!C Calculate the multi-body contribution to energy.
+! ecorr=ecorr+ekont*ees
+!C Calculate multi-body contributions to the gradient.
+ coeffpees0pij=coeffp*ees0pij
+ coeffmees0mij=coeffm*ees0mij
+ coeffpees0pkl=coeffp*ees0pkl
+ coeffmees0mkl=coeffm*ees0mkl
+ do ll=1,3
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i) &
+ -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+&
+ coeffmees0mkl*gacontm_hb1(ll,jj,i))
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j) &
+ -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb2(ll,jj,i))
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k) &
+ -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb1(ll,kk,k))
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l) &
+ -ekont*(coeffpees0pij*gacontp_hb2(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb2(ll,kk,k))
+ gradlongij=ees*ekl*gacont_hbr(ll,jj,i)- &
+ ekont*(coeffpees0pkl*gacontp_hb3(ll,jj,i)+ &
+ coeffmees0mkl*gacontm_hb3(ll,jj,i))
+ gradcorr3_nucl(ll,j)=gradcorr3_nucl(ll,j)+gradlongij
+ gradcorr3_nucl(ll,i)=gradcorr3_nucl(ll,i)-gradlongij
+ gradlongkl=ees*eij*gacont_hbr(ll,kk,k)- &
+ ekont*(coeffpees0pij*gacontp_hb3(ll,kk,k)+ &
+ coeffmees0mij*gacontm_hb3(ll,kk,k))
+ gradcorr3_nucl(ll,l)=gradcorr3_nucl(ll,l)+gradlongkl
+ gradcorr3_nucl(ll,k)=gradcorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,i)=gradxorr3_nucl(ll,i)-gradlongij
+ gradxorr3_nucl(ll,j)=gradxorr3_nucl(ll,j)+gradlongij
+ gradxorr3_nucl(ll,k)=gradxorr3_nucl(ll,k)-gradlongkl
+ gradxorr3_nucl(ll,l)=gradxorr3_nucl(ll,l)+gradlongkl
+ enddo
+ ehbcorr3_nucl=ekont*ees
+ return
+ end function ehbcorr3_nucl
+#ifdef MPI
+ subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+ num_kont=num_cont_hb(atom)
+ do i=1,num_kont
+ do k=1,8
+ do j=1,3
+ buffer(i,indx+(k-1)*3+j)=zapas2(j,i,atom,k)
+ enddo ! j
+ enddo ! k
+ buffer(i,indx+25)=facont_hb(i,atom)
+ buffer(i,indx+26)=ees0p(i,atom)
+ buffer(i,indx+27)=ees0m(i,atom)
+ buffer(i,indx+28)=d_cont(i,atom)
+ buffer(i,indx+29)=dfloat(jcont_hb(i,atom))
+ enddo ! i
+ buffer(1,indx+30)=dfloat(num_kont)
+ return
+ end subroutine pack_buffer
+!c------------------------------------------------------------------------------
+ subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer)
+ integer dimen1,dimen2,atom,indx,numcont,i,ii,k,j,num_kont,num_kont_old
+ real(kind=8):: buffer(dimen1,dimen2)
+! double precision zapas
+! common /contacts_hb/ zapas(3,maxconts,maxres,8),
+! & facont_hb(maxconts,maxres),ees0p(maxconts,maxres),
+! & ees0m(maxconts,maxres),d_cont(maxconts,maxres),
+! & num_cont_hb(maxres),jcont_hb(maxconts,maxres)
+ num_kont=buffer(1,indx+30)
+ num_kont_old=num_cont_hb(atom)
+ num_cont_hb(atom)=num_kont+num_kont_old
+ do i=1,num_kont
+ ii=i+num_kont_old
+ do k=1,8
+ do j=1,3
+ zapas2(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j)
+ enddo ! j
+ enddo ! k
+ facont_hb(ii,atom)=buffer(i,indx+25)
+ ees0p(ii,atom)=buffer(i,indx+26)
+ ees0m(ii,atom)=buffer(i,indx+27)
+ d_cont(i,atom)=buffer(i,indx+28)
+ jcont_hb(ii,atom)=buffer(i,indx+29)
+ enddo ! i
+ return
+ end subroutine unpack_buffer
+!c------------------------------------------------------------------------------
+#endif
!----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
projects / unres4.git / blobdiff