From 3c2962eb07d557194f5368c06ed6a574a2f16e3b Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Fri, 11 Aug 2017 13:48:50 +0200 Subject: [PATCH] energi calculation working (no weights!) --- source/unres/control.F90 | 2 +- source/unres/data/energy_data.f90 | 1 + source/unres/energy.f90 | 441 +++++++++++++++++++++++++++++++++++-- 3 files changed, 430 insertions(+), 14 deletions(-) diff --git a/source/unres/control.F90 b/source/unres/control.F90 index f8f8e9a..d70d144 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -509,7 +509,7 @@ my_ele_inds_vdw_nucl,my_ele_inde_vdw_nucl,ind_eleint_vdw_nucl,& ind_eleint_vdw_old_nucl,nscp_int_tot_nucl,my_scp_inds_nucl,& my_scp_inde_nucl,ind_scpint_nucl,ind_scpint_old_nucl - integer,dimension(5) :: nct_molec,nnt_molec +! integer,dimension(5) :: nct_molec,nnt_molec !el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1) !el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1) diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index 97b773c..bbc2ed2 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -104,6 +104,7 @@ real(kind=8),dimension(:,:),allocatable :: ael6_nucl,& ael3_nucl,ael32_nucl,ael63_nucl integer :: expon,expon2, nnt,nct,itypro + integer,dimension(5) :: nnt_molec,nct_molec integer,dimension(:,:),allocatable :: istart,iend !(maxres,maxint_gr) integer,dimension(:),allocatable :: nint_gr,itel,& ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr !(maxres) diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index b20eb96..881497c 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -128,7 +128,7 @@ !-----------------------------NUCLEIC GRADIENT real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl, & gvdwpsb1,gelpp,gvdwpsb,gelsbc,gelsbx,gvdwsbx,gvdwsbc,gsbloc,& - gsblocx + 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) @@ -185,7 +185,8 @@ 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) !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- @@ -548,14 +549,15 @@ 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 @@ -3065,6 +3067,7 @@ ! ! 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) @@ -3196,10 +3199,10 @@ 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) @@ -19597,6 +19600,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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)) @@ -20412,7 +20419,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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,& @@ -20428,6 +20435,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! 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 @@ -20572,28 +20580,28 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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/ @@ -20716,7 +20724,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 @@ -20729,6 +20738,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !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.' @@ -21075,6 +21086,410 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 !---------------------------------------------------------------------------- !----------------------------------------------------------------------------- -- 1.7.9.5