1c1 < subroutine etotal(energia,fact) --- > subroutine etotal(energia) 4,5d3 < include 'DIMENSIONS.ZSCOPT' < 8d5 < #endif 12,18d8 < < include 'COMMON.IOUNITS' < double precision energia(0:max_ene),energia1(0:max_ene+1) < #ifdef MPL < include 'COMMON.INFO' < external d_vadd < integer ready 19a10,17 > #ifdef MPI > include "mpif.h" > double precision weights_(n_ene) > #endif > include 'COMMON.SETUP' > include 'COMMON.IOUNITS' > double precision energia(0:n_ene) > include 'COMMON.LOCAL' 25,28c23,98 < double precision fact(6) < cd write(iout, '(a,i2)')'Calling etotal ipot=',ipot < cd print *,'nnt=',nnt,' nct=',nct < C --- > include 'COMMON.VAR' > include 'COMMON.MD' > include 'COMMON.CONTROL' > include 'COMMON.TIME1' > #ifdef MPI > c print*,"ETOTAL Processor",fg_rank," absolute rank",myrank, > c & " nfgtasks",nfgtasks > if (nfgtasks.gt.1) then > time00=MPI_Wtime() > C FG slaves call the following matching MPI_Bcast in ERGASTULUM > if (fg_rank.eq.0) then > call MPI_Bcast(0,1,MPI_INTEGER,king,FG_COMM,IERROR) > c print *,"Processor",myrank," BROADCAST iorder" > C FG master sets up the WEIGHTS_ array which will be broadcast to the > C FG slaves as WEIGHTS array. > weights_(1)=wsc > weights_(2)=wscp > weights_(3)=welec > weights_(4)=wcorr > weights_(5)=wcorr5 > weights_(6)=wcorr6 > weights_(7)=wel_loc > weights_(8)=wturn3 > weights_(9)=wturn4 > weights_(10)=wturn6 > weights_(11)=wang > weights_(12)=wscloc > weights_(13)=wtor > weights_(14)=wtor_d > weights_(15)=wstrain > weights_(16)=wvdwpp > weights_(17)=wbond > weights_(18)=scal14 > weights_(21)=wsccor > C FG Master broadcasts the WEIGHTS_ array > call MPI_Bcast(weights_(1),n_ene, > & MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) > else > C FG slaves receive the WEIGHTS array > call MPI_Bcast(weights(1),n_ene, > & MPI_DOUBLE_PRECISION,king,FG_COMM,IERROR) > wsc=weights(1) > wscp=weights(2) > welec=weights(3) > wcorr=weights(4) > wcorr5=weights(5) > wcorr6=weights(6) > wel_loc=weights(7) > wturn3=weights(8) > wturn4=weights(9) > wturn6=weights(10) > wang=weights(11) > wscloc=weights(12) > wtor=weights(13) > wtor_d=weights(14) > wstrain=weights(15) > wvdwpp=weights(16) > wbond=weights(17) > scal14=weights(18) > wsccor=weights(21) > endif > time_Bcast=time_Bcast+MPI_Wtime()-time00 > time_Bcastw=time_Bcastw+MPI_Wtime()-time00 > c call chainbuild_cart > endif > c print *,'Processor',myrank,' calling etotal ipot=',ipot > c print *,'Processor',myrank,' nnt=',nnt,' nct=',nct > #else > c if (modecalc.eq.12.or.modecalc.eq.14) then > c call int_from_cart1(.false.) > c endif > #endif > #ifdef TIMING > time00=MPI_Wtime() > #endif > C 31c101 < goto (101,102,103,104,105) ipot --- > goto (101,102,103,104,105,106) ipot 33c103 < 101 call elj(evdw,evdw_t) --- > 101 call elj(evdw) 35c105 < goto 106 --- > goto 107 37,38c107,108 < 102 call eljk(evdw,evdw_t) < goto 106 --- > 102 call eljk(evdw) > goto 107 40,41c110,111 < 103 call ebp(evdw,evdw_t) < goto 106 --- > 103 call ebp(evdw) > goto 107 43,44c113,114 < 104 call egb(evdw,evdw_t) < goto 106 --- > 104 call egb(evdw) > goto 107 46c116,119 < 105 call egbv(evdw,evdw_t) --- > 105 call egbv(evdw) > goto 107 > C Soft-sphere potential > 106 call e_softsphere(evdw) 50c123,158 < 106 call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) --- > 107 continue > c print *,"Processor",myrank," computed USCSC" > #ifdef TIMING > time01=MPI_Wtime() > #endif > call vec_and_deriv > #ifdef TIMING > time_vec=time_vec+MPI_Wtime()-time01 > #endif > c print *,"Processor",myrank," left VEC_AND_DERIV" > if (ipot.lt.6) then > #ifdef SPLITELE > if (welec.gt.0d0.or.wvdwpp.gt.0d0.or.wel_loc.gt.0d0.or. > & wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 > & .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 > & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then > #else > if (welec.gt.0d0.or.wel_loc.gt.0d0.or. > & wturn3.gt.0d0.or.wturn4.gt.0d0 .or. wcorr.gt.0.0d0 > & .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 > & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then > #endif > call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) > else > ees=0.0d0 > evdw1=0.0d0 > eel_loc=0.0d0 > eello_turn3=0.0d0 > eello_turn4=0.0d0 > endif > else > c write (iout,*) "Soft-spheer ELEC potential" > call eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, > & eello_turn4) > endif > c print *,"Processor",myrank," computed UELEC" 55c163,173 < call escp(evdw2,evdw2_14) --- > if (ipot.lt.6) then > if(wscp.gt.0d0) then > call escp(evdw2,evdw2_14) > else > evdw2=0 > evdw2_14=0 > endif > else > c write (iout,*) "Soft-sphere SCP potential" > call escp_soft_sphere(evdw2,evdw2_14) > endif 60d177 < c write (iout,*) "estr",estr 70,71c187,192 < call ebend(ebe) < cd print *,'Bend energy finished.' --- > if (wang.gt.0d0) then > call ebend(ebe) > else > ebe=0 > endif > c print *,"Processor",myrank," computed UB" 76c197 < cd print *,'SCLOC energy finished.' --- > c print *,"Processor",myrank," computed USC" 81c202,208 < call etor(etors,edihcnstr,fact(1)) --- > if (wtor.gt.0) then > call etor(etors,edihcnstr) > else > etors=0 > edihcnstr=0 > endif > c print *,"Processor",myrank," computed Utor" 85c212,217 < call etor_d(etors_d,fact(2)) --- > if (wtor_d.gt.0) then > call etor_d(etors_d) > else > etors_d=0 > endif > c print *,"Processor",myrank," computed Utord" 89c221,226 < call eback_sc_corr(esccor) --- > if (wsccor.gt.0.0d0) then > call eback_sc_corr(esccor) > else > esccor=0.0d0 > endif > c print *,"Processor",myrank," computed Usccorr" 95,97c232,233 < if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 < & .or. wturn6.gt.0.0d0) then < c print *,"calling multibody_eello" --- > if ((wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 > & .or. wturn6.gt.0.0d0) .and. ipot.lt.6) then 99,100c235,241 < c write (*,*) 'n_corr=',n_corr,' n_corr1=',n_corr1 < c print *,ecorr,ecorr5,ecorr6,eturn6 --- > cd write(2,*)'multibody_eello n_corr=',n_corr,' n_corr1=',n_corr1, > cd &" ecorr",ecorr," ecorr5",ecorr5," ecorr6",ecorr6," eturn6",eturn6 > else > ecorr=0.0d0 > ecorr5=0.0d0 > ecorr6=0.0d0 > eturn6=0.0d0 102c243 < if (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) then --- > if ((wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) .and. ipot.lt.6) then 103a245 > cd write (iout,*) "multibody_hb ecorr",ecorr 105,123c247,259 < c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t < #ifdef SPLITELE < etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees < & +wvdwpp*evdw1 < & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc < & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 < & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 < & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 < & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d < & +wbond*estr+wsccor*fact(1)*esccor < #else < etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 < & +welec*fact(1)*(ees+evdw1) < & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc < & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 < & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 < & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 < & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d < & +wbond*estr+wsccor*fact(1)*esccor --- > c print *,"Processor",myrank," computed Ucorr" > C > C If performing constraint dynamics, call the constraint energy > C after the equilibration time > if(usampl.and.totT.gt.eq_time) then > call EconstrQ > call Econstr_back > else > Uconst=0.0d0 > Uconst_back=0.0d0 > endif > #ifdef TIMING > time_enecalc=time_enecalc+MPI_Wtime()-time00 125c261,267 < energia(0)=etot --- > c print *,"Processor",myrank," computed Uconstr" > #ifdef TIMING > time00=MPI_Wtime() > #endif > c > C Sum the energies > C 129c271 < energia(17)=evdw2_14 --- > energia(18)=evdw2_14 132c274 < energia(17)=0.0d0 --- > energia(18)=0.0d0 153,156c295,402 < energia(18)=estr < energia(19)=esccor < energia(20)=edihcnstr < energia(21)=evdw_t --- > energia(19)=edihcnstr > energia(17)=estr > energia(20)=Uconst+Uconst_back > energia(21)=esccor > c print *," Processor",myrank," calls SUM_ENERGY" > call sum_energy(energia,.true.) > c print *," Processor",myrank," left SUM_ENERGY" > #ifdef TIMING > time_sumene=time_sumene+MPI_Wtime()-time00 > #endif > return > end > c------------------------------------------------------------------------------- > subroutine sum_energy(energia,reduce) > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > #ifndef ISNAN > external proc_proc > #ifdef WINPGI > cMS$ATTRIBUTES C :: proc_proc > #endif > #endif > #ifdef MPI > include "mpif.h" > #endif > include 'COMMON.SETUP' > include 'COMMON.IOUNITS' > double precision energia(0:n_ene),enebuff(0:n_ene+1) > include 'COMMON.FFIELD' > include 'COMMON.DERIV' > include 'COMMON.INTERACT' > include 'COMMON.SBRIDGE' > include 'COMMON.CHAIN' > include 'COMMON.VAR' > include 'COMMON.CONTROL' > include 'COMMON.TIME1' > logical reduce > #ifdef MPI > if (nfgtasks.gt.1 .and. reduce) then > #ifdef DEBUG > write (iout,*) "energies before REDUCE" > call enerprint(energia) > call flush(iout) > #endif > do i=0,n_ene > enebuff(i)=energia(i) > enddo > time00=MPI_Wtime() > call MPI_Barrier(FG_COMM,IERR) > time_barrier_e=time_barrier_e+MPI_Wtime()-time00 > time00=MPI_Wtime() > call MPI_Reduce(enebuff(0),energia(0),n_ene+1, > & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) > #ifdef DEBUG > write (iout,*) "energies after REDUCE" > call enerprint(energia) > call flush(iout) > #endif > time_Reduce=time_Reduce+MPI_Wtime()-time00 > endif > if (fg_rank.eq.0) then > #endif > evdw=energia(1) > #ifdef SCP14 > evdw2=energia(2)+energia(18) > evdw2_14=energia(18) > #else > evdw2=energia(2) > #endif > #ifdef SPLITELE > ees=energia(3) > evdw1=energia(16) > #else > ees=energia(3) > evdw1=0.0d0 > #endif > ecorr=energia(4) > ecorr5=energia(5) > ecorr6=energia(6) > eel_loc=energia(7) > eello_turn3=energia(8) > eello_turn4=energia(9) > eturn6=energia(10) > ebe=energia(11) > escloc=energia(12) > etors=energia(13) > etors_d=energia(14) > ehpb=energia(15) > edihcnstr=energia(19) > estr=energia(17) > Uconst=energia(20) > esccor=energia(21) > #ifdef SPLITELE > etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 > & +wang*ebe+wtor*etors+wscloc*escloc > & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 > & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 > & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d > & +wbond*estr+Uconst+wsccor*esccor > #else > etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) > & +wang*ebe+wtor*etors+wscloc*escloc > & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 > & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 > & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d > & +wbond*estr+Uconst+wsccor*esccor > #endif > energia(0)=etot 173,174c419,464 < #ifdef MPL < c endif --- > #ifdef MPI > endif > #endif > return > end > c------------------------------------------------------------------------------- > subroutine sum_gradient > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > #ifndef ISNAN > external proc_proc > #ifdef WINPGI > cMS$ATTRIBUTES C :: proc_proc > #endif > #endif > #ifdef MPI > include 'mpif.h' > double precision gradbufc(3,maxres),gradbufx(3,maxres), > & glocbuf(4*maxres),gradbufc_sum(3,maxres) > #endif > include 'COMMON.SETUP' > include 'COMMON.IOUNITS' > include 'COMMON.FFIELD' > include 'COMMON.DERIV' > include 'COMMON.INTERACT' > include 'COMMON.SBRIDGE' > include 'COMMON.CHAIN' > include 'COMMON.VAR' > include 'COMMON.CONTROL' > include 'COMMON.TIME1' > include 'COMMON.MAXGRAD' > #ifdef TIMING > time01=MPI_Wtime() > #endif > #ifdef DEBUG > write (iout,*) "sum_gradient gvdwc, gvdwx" > do i=1,nres > write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') > & i,(gvdwx(j,i),j=1,3),(gvdwc(j,i),j=1,3) > enddo > call flush(iout) > #endif > #ifdef MPI > C FG slaves call the following matching MPI_Bcast in ERGASTULUM > if (nfgtasks.gt.1 .and. fg_rank.eq.0) > & call MPI_Bcast(1,1,MPI_INTEGER,king,FG_COMM,IERROR) 176d465 < if (calc_grad) then 178c467,468 < C Sum up the components of the Cartesian gradient. --- > C 9/29/08 AL Transform parts of gradients in site coordinates to the gradient > C in virtual-bond-vector coordinates 179a470,488 > #ifdef DEBUG > c write (iout,*) "gel_loc gel_loc_long and gel_loc_loc" > c do i=1,nres-1 > c write (iout,'(i5,3f10.5,2x,3f10.5,2x,f10.5)') > c & i,(gel_loc(j,i),j=1,3),(gel_loc_long(j,i),j=1,3),gel_loc_loc(i) > c enddo > c write (iout,*) "gel_loc_tur3 gel_loc_turn4" > c do i=1,nres-1 > c write (iout,'(i5,3f10.5,2x,f10.5)') > c & i,(gcorr4_turn(j,i),j=1,3),gel_loc_turn4(i) > c enddo > write (iout,*) "gradcorr5 gradcorr5_long gradcorr5_loc" > do i=1,nres > write (iout,'(i3,3f10.5,5x,3f10.5,5x,f10.5)') > & i,(gradcorr5(j,i),j=1,3),(gradcorr5_long(j,i),j=1,3), > & g_corr5_loc(i) > enddo > call flush(iout) > #endif 183,198c492,500 < gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ < & welec*fact(1)*gelc(j,i)+wvdwpp*gvdwpp(j,i)+ < & wbond*gradb(j,i)+ < & wstrain*ghpbc(j,i)+ < & wcorr*fact(3)*gradcorr(j,i)+ < & wel_loc*fact(2)*gel_loc(j,i)+ < & wturn3*fact(2)*gcorr3_turn(j,i)+ < & wturn4*fact(3)*gcorr4_turn(j,i)+ < & wcorr5*fact(4)*gradcorr5(j,i)+ < & wcorr6*fact(5)*gradcorr6(j,i)+ < & wturn6*fact(5)*gcorr6_turn(j,i)+ < & wsccor*fact(2)*gsccorc(j,i) < gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ < & wbond*gradbx(j,i)+ < & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ < & wsccor*fact(2)*gsccorx(j,i) --- > gradbufc(j,i)=wsc*gvdwc(j,i)+ > & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ > & welec*gelc_long(j,i)+wvdwpp*gvdwpp(j,i)+ > & wel_loc*gel_loc_long(j,i)+ > & wcorr*gradcorr_long(j,i)+ > & wcorr5*gradcorr5_long(j,i)+ > & wcorr6*gradcorr6_long(j,i)+ > & wturn6*gcorr6_turn_long(j,i)+ > & wstrain*ghpbc(j,i) 199a502 > enddo 203,204c506,508 < gradc(j,i,icg)=wsc*gvdwc(j,i)+wscp*gvdwc_scp(j,i)+ < & welec*fact(1)*gelc(j,i)+wstrain*ghpbc(j,i)+ --- > gradbufc(j,i)=wsc*gvdwc(j,i)+ > & wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+ > & welec*gelc_long(j,i)+ 206,213c510,670 < & wcorr*fact(3)*gradcorr(j,i)+ < & wel_loc*fact(2)*gel_loc(j,i)+ < & wturn3*fact(2)*gcorr3_turn(j,i)+ < & wturn4*fact(3)*gcorr4_turn(j,i)+ < & wcorr5*fact(4)*gradcorr5(j,i)+ < & wcorr6*fact(5)*gradcorr6(j,i)+ < & wturn6*fact(5)*gcorr6_turn(j,i)+ < & wsccor*fact(2)*gsccorc(j,i) --- > & wel_loc*gel_loc_long(j,i)+ > & wcorr*gradcorr_long(j,i)+ > & wcorr5*gradcorr5_long(j,i)+ > & wcorr6*gradcorr6_long(j,i)+ > & wturn6*gcorr6_turn_long(j,i)+ > & wstrain*ghpbc(j,i) > enddo > enddo > #endif > #ifdef MPI > if (nfgtasks.gt.1) then > time00=MPI_Wtime() > #ifdef DEBUG > write (iout,*) "gradbufc before allreduce" > do i=1,nres > write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) > enddo > call flush(iout) > #endif > do i=1,nres > do j=1,3 > gradbufc_sum(j,i)=gradbufc(j,i) > enddo > enddo > c call MPI_AllReduce(gradbufc(1,1),gradbufc_sum(1,1),3*nres, > c & MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR) > c time_reduce=time_reduce+MPI_Wtime()-time00 > #ifdef DEBUG > c write (iout,*) "gradbufc_sum after allreduce" > c do i=1,nres > c write (iout,'(i3,3f10.5)') i,(gradbufc_sum(j,i),j=1,3) > c enddo > c call flush(iout) > #endif > #ifdef TIMING > c time_allreduce=time_allreduce+MPI_Wtime()-time00 > #endif > do i=nnt,nres > do k=1,3 > gradbufc(k,i)=0.0d0 > enddo > enddo > #ifdef DEBUG > write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end > write (iout,*) (i," jgrad_start",jgrad_start(i), > & " jgrad_end ",jgrad_end(i), > & i=igrad_start,igrad_end) > #endif > c > c Obsolete and inefficient code; we can make the effort O(n) and, therefore, > c do not parallelize this part. > c > c do i=igrad_start,igrad_end > c do j=jgrad_start(i),jgrad_end(i) > c do k=1,3 > c gradbufc(k,i)=gradbufc(k,i)+gradbufc_sum(k,j) > c enddo > c enddo > c enddo > do j=1,3 > gradbufc(j,nres-1)=gradbufc_sum(j,nres) > enddo > do i=nres-2,nnt,-1 > do j=1,3 > gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) > enddo > enddo > #ifdef DEBUG > write (iout,*) "gradbufc after summing" > do i=1,nres > write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) > enddo > call flush(iout) > #endif > else > #endif > #ifdef DEBUG > write (iout,*) "gradbufc" > do i=1,nres > write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) > enddo > call flush(iout) > #endif > do i=1,nres > do j=1,3 > gradbufc_sum(j,i)=gradbufc(j,i) > gradbufc(j,i)=0.0d0 > enddo > enddo > do j=1,3 > gradbufc(j,nres-1)=gradbufc_sum(j,nres) > enddo > do i=nres-2,nnt,-1 > do j=1,3 > gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1) > enddo > enddo > c do i=nnt,nres-1 > c do k=1,3 > c gradbufc(k,i)=0.0d0 > c enddo > c do j=i+1,nres > c do k=1,3 > c gradbufc(k,i)=gradbufc(k,i)+gradbufc(k,j) > c enddo > c enddo > c enddo > #ifdef DEBUG > write (iout,*) "gradbufc after summing" > do i=1,nres > write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3) > enddo > call flush(iout) > #endif > #ifdef MPI > endif > #endif > do k=1,3 > gradbufc(k,nres)=0.0d0 > enddo > do i=1,nct > do j=1,3 > #ifdef SPLITELE > gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ > & wel_loc*gel_loc(j,i)+ > & 0.5d0*(wscp*gvdwc_scpp(j,i)+ > & welec*gelc_long(j,i)+wvdwpp*gvdwpp(j,i)+ > & wel_loc*gel_loc_long(j,i)+ > & wcorr*gradcorr_long(j,i)+ > & wcorr5*gradcorr5_long(j,i)+ > & wcorr6*gradcorr6_long(j,i)+ > & wturn6*gcorr6_turn_long(j,i))+ > & wbond*gradb(j,i)+ > & wcorr*gradcorr(j,i)+ > & wturn3*gcorr3_turn(j,i)+ > & wturn4*gcorr4_turn(j,i)+ > & wcorr5*gradcorr5(j,i)+ > & wcorr6*gradcorr6(j,i)+ > & wturn6*gcorr6_turn(j,i)+ > & wsccor*gsccorc(j,i) > & +wscloc*gscloc(j,i) > #else > gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ > & wel_loc*gel_loc(j,i)+ > & 0.5d0*(wscp*gvdwc_scpp(j,i)+ > & welec*gelc_long(j,i) > & wel_loc*gel_loc_long(j,i)+ > & wcorr*gcorr_long(j,i)+ > & wcorr5*gradcorr5_long(j,i)+ > & wcorr6*gradcorr6_long(j,i)+ > & wturn6*gcorr6_turn_long(j,i))+ > & wbond*gradb(j,i)+ > & wcorr*gradcorr(j,i)+ > & wturn3*gcorr3_turn(j,i)+ > & wturn4*gcorr4_turn(j,i)+ > & wcorr5*gradcorr5(j,i)+ > & wcorr6*gradcorr6(j,i)+ > & wturn6*gcorr6_turn(j,i)+ > & wsccor*gsccorc(j,i) > & +wscloc*gscloc(j,i) > #endif 217c674,675 < & wsccor*fact(1)*gsccorx(j,i) --- > & wsccor*gsccorx(j,i) > & +wscloc*gsclocx(j,i) 219c677,681 < #endif --- > enddo > #ifdef DEBUG > write (iout,*) "gloc before adding corr" > do i=1,4*nres > write (iout,*) i,gloc(i,icg) 221,222c683 < < --- > #endif 224,231c685,697 < gloc(i,icg)=gloc(i,icg)+wcorr*fact(3)*gcorr_loc(i) < & +wcorr5*fact(4)*g_corr5_loc(i) < & +wcorr6*fact(5)*g_corr6_loc(i) < & +wturn4*fact(3)*gel_loc_turn4(i) < & +wturn3*fact(2)*gel_loc_turn3(i) < & +wturn6*fact(5)*gel_loc_turn6(i) < & +wel_loc*fact(2)*gel_loc_loc(i) < & +wsccor*fact(1)*gsccor_loc(i) --- > gloc(i,icg)=gloc(i,icg)+wcorr*gcorr_loc(i) > & +wcorr5*g_corr5_loc(i) > & +wcorr6*g_corr6_loc(i) > & +wturn4*gel_loc_turn4(i) > & +wturn3*gel_loc_turn3(i) > & +wturn6*gel_loc_turn6(i) > & +wel_loc*gel_loc_loc(i) > & +wsccor*gsccor_loc(i) > enddo > #ifdef DEBUG > write (iout,*) "gloc after adding corr" > do i=1,4*nres > write (iout,*) i,gloc(i,icg) 232a699,727 > #endif > #ifdef MPI > if (nfgtasks.gt.1) then > do j=1,3 > do i=1,nres > gradbufc(j,i)=gradc(j,i,icg) > gradbufx(j,i)=gradx(j,i,icg) > enddo > enddo > do i=1,4*nres > glocbuf(i)=gloc(i,icg) > enddo > time00=MPI_Wtime() > call MPI_Barrier(FG_COMM,IERR) > time_barrier_g=time_barrier_g+MPI_Wtime()-time00 > time00=MPI_Wtime() > call MPI_Reduce(gradbufc(1,1),gradc(1,1,icg),3*nres, > & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) > call MPI_Reduce(gradbufx(1,1),gradx(1,1,icg),3*nres, > & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) > call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres, > & MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR) > time_reduce=time_reduce+MPI_Wtime()-time00 > #ifdef DEBUG > write (iout,*) "gloc after reduce" > do i=1,4*nres > write (iout,*) i,gloc(i,icg) > enddo > #endif 233a729,897 > #endif > if (gnorm_check) then > c > c Compute the maximum elements of the gradient > c > gvdwc_max=0.0d0 > gvdwc_scp_max=0.0d0 > gelc_max=0.0d0 > gvdwpp_max=0.0d0 > gradb_max=0.0d0 > ghpbc_max=0.0d0 > gradcorr_max=0.0d0 > gel_loc_max=0.0d0 > gcorr3_turn_max=0.0d0 > gcorr4_turn_max=0.0d0 > gradcorr5_max=0.0d0 > gradcorr6_max=0.0d0 > gcorr6_turn_max=0.0d0 > gsccorc_max=0.0d0 > gscloc_max=0.0d0 > gvdwx_max=0.0d0 > gradx_scp_max=0.0d0 > ghpbx_max=0.0d0 > gradxorr_max=0.0d0 > gsccorx_max=0.0d0 > gsclocx_max=0.0d0 > do i=1,nct > gvdwc_norm=dsqrt(scalar(gvdwc(1,i),gvdwc(1,i))) > if (gvdwc_norm.gt.gvdwc_max) gvdwc_max=gvdwc_norm > gvdwc_scp_norm=dsqrt(scalar(gvdwc_scp(1,i),gvdwc_scp(1,i))) > if (gvdwc_scp_norm.gt.gvdwc_scp_max) > & gvdwc_scp_max=gvdwc_scp_norm > gelc_norm=dsqrt(scalar(gelc(1,i),gelc(1,i))) > if (gelc_norm.gt.gelc_max) gelc_max=gelc_norm > gvdwpp_norm=dsqrt(scalar(gvdwpp(1,i),gvdwpp(1,i))) > if (gvdwpp_norm.gt.gvdwpp_max) gvdwpp_max=gvdwpp_norm > gradb_norm=dsqrt(scalar(gradb(1,i),gradb(1,i))) > if (gradb_norm.gt.gradb_max) gradb_max=gradb_norm > ghpbc_norm=dsqrt(scalar(ghpbc(1,i),ghpbc(1,i))) > if (ghpbc_norm.gt.ghpbc_max) ghpbc_max=ghpbc_norm > gradcorr_norm=dsqrt(scalar(gradcorr(1,i),gradcorr(1,i))) > if (gradcorr_norm.gt.gradcorr_max) gradcorr_max=gradcorr_norm > gel_loc_norm=dsqrt(scalar(gel_loc(1,i),gel_loc(1,i))) > if (gel_loc_norm.gt.gel_loc_max) gel_loc_max=gel_loc_norm > gcorr3_turn_norm=dsqrt(scalar(gcorr3_turn(1,i), > & gcorr3_turn(1,i))) > if (gcorr3_turn_norm.gt.gcorr3_turn_max) > & gcorr3_turn_max=gcorr3_turn_norm > gcorr4_turn_norm=dsqrt(scalar(gcorr4_turn(1,i), > & gcorr4_turn(1,i))) > if (gcorr4_turn_norm.gt.gcorr4_turn_max) > & gcorr4_turn_max=gcorr4_turn_norm > gradcorr5_norm=dsqrt(scalar(gradcorr5(1,i),gradcorr5(1,i))) > if (gradcorr5_norm.gt.gradcorr5_max) > & gradcorr5_max=gradcorr5_norm > gradcorr6_norm=dsqrt(scalar(gradcorr6(1,i),gradcorr6(1,i))) > if (gradcorr6_norm.gt.gradcorr6_max) gcorr6_max=gradcorr6_norm > gcorr6_turn_norm=dsqrt(scalar(gcorr6_turn(1,i), > & gcorr6_turn(1,i))) > if (gcorr6_turn_norm.gt.gcorr6_turn_max) > & gcorr6_turn_max=gcorr6_turn_norm > gsccorr_norm=dsqrt(scalar(gsccorc(1,i),gsccorc(1,i))) > if (gsccorr_norm.gt.gsccorr_max) gsccorr_max=gsccorr_norm > gscloc_norm=dsqrt(scalar(gscloc(1,i),gscloc(1,i))) > if (gscloc_norm.gt.gscloc_max) gscloc_max=gscloc_norm > gvdwx_norm=dsqrt(scalar(gvdwx(1,i),gvdwx(1,i))) > if (gvdwx_norm.gt.gvdwx_max) gvdwx_max=gvdwx_norm > gradx_scp_norm=dsqrt(scalar(gradx_scp(1,i),gradx_scp(1,i))) > if (gradx_scp_norm.gt.gradx_scp_max) > & gradx_scp_max=gradx_scp_norm > ghpbx_norm=dsqrt(scalar(ghpbx(1,i),ghpbx(1,i))) > if (ghpbx_norm.gt.ghpbx_max) ghpbx_max=ghpbx_norm > gradxorr_norm=dsqrt(scalar(gradxorr(1,i),gradxorr(1,i))) > if (gradxorr_norm.gt.gradxorr_max) gradxorr_max=gradxorr_norm > gsccorrx_norm=dsqrt(scalar(gsccorx(1,i),gsccorx(1,i))) > if (gsccorrx_norm.gt.gsccorrx_max) gsccorrx_max=gsccorrx_norm > gsclocx_norm=dsqrt(scalar(gsclocx(1,i),gsclocx(1,i))) > if (gsclocx_norm.gt.gsclocx_max) gsclocx_max=gsclocx_norm > enddo > if (gradout) then > #ifdef AIX > open(istat,file=statname,position="append") > #else > open(istat,file=statname,access="append") > #endif > write (istat,'(1h#,21f10.2)') gvdwc_max,gvdwc_scp_max, > & gelc_max,gvdwpp_max,gradb_max,ghpbc_max, > & gradcorr_max,gel_loc_max,gcorr3_turn_max,gcorr4_turn_max, > & gradcorr5_max,gradcorr6_max,gcorr6_turn_max,gsccorc_max, > & gscloc_max,gvdwx_max,gradx_scp_max,ghpbx_max,gradxorr_max, > & gsccorx_max,gsclocx_max > close(istat) > if (gvdwc_max.gt.1.0d4) then > write (iout,*) "gvdwc gvdwx gradb gradbx" > do i=nnt,nct > write(iout,'(i5,4(3f10.2,5x))') i,(gvdwc(j,i),gvdwx(j,i), > & gradb(j,i),gradbx(j,i),j=1,3) > enddo > call pdbout(0.0d0,'cipiszcze',iout) > call flush(iout) > endif > endif > endif > #ifdef DEBUG > write (iout,*) "gradc gradx gloc" > do i=1,nres > write (iout,'(i5,3f10.5,5x,3f10.5,5x,f10.5)') > & i,(gradc(j,i,icg),j=1,3),(gradx(j,i,icg),j=1,3),gloc(i,icg) > enddo > #endif > #ifdef TIMING > time_sumgradient=time_sumgradient+MPI_Wtime()-time01 > #endif > return > end > c------------------------------------------------------------------------------- > subroutine rescale_weights(t_bath) > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'COMMON.IOUNITS' > include 'COMMON.FFIELD' > include 'COMMON.SBRIDGE' > double precision kfac /2.4d0/ > double precision x,x2,x3,x4,x5,licznik /1.12692801104297249644/ > c facT=temp0/t_bath > c facT=2*temp0/(t_bath+temp0) > if (rescale_mode.eq.0) then > facT=1.0d0 > facT2=1.0d0 > facT3=1.0d0 > facT4=1.0d0 > facT5=1.0d0 > else if (rescale_mode.eq.1) then > facT=kfac/(kfac-1.0d0+t_bath/temp0) > facT2=kfac**2/(kfac**2-1.0d0+(t_bath/temp0)**2) > facT3=kfac**3/(kfac**3-1.0d0+(t_bath/temp0)**3) > facT4=kfac**4/(kfac**4-1.0d0+(t_bath/temp0)**4) > facT5=kfac**5/(kfac**5-1.0d0+(t_bath/temp0)**5) > else if (rescale_mode.eq.2) then > x=t_bath/temp0 > x2=x*x > x3=x2*x > x4=x3*x > x5=x4*x > facT=licznik/dlog(dexp(x)+dexp(-x)) > facT2=licznik/dlog(dexp(x2)+dexp(-x2)) > facT3=licznik/dlog(dexp(x3)+dexp(-x3)) > facT4=licznik/dlog(dexp(x4)+dexp(-x4)) > facT5=licznik/dlog(dexp(x5)+dexp(-x5)) > else > write (iout,*) "Wrong RESCALE_MODE",rescale_mode > write (*,*) "Wrong RESCALE_MODE",rescale_mode > #ifdef MPI > call MPI_Finalize(MPI_COMM_WORLD,IERROR) > #endif > stop 555 > endif > welec=weights(3)*fact > wcorr=weights(4)*fact3 > wcorr5=weights(5)*fact4 > wcorr6=weights(6)*fact5 > wel_loc=weights(7)*fact2 > wturn3=weights(8)*fact2 > wturn4=weights(9)*fact3 > wturn6=weights(10)*fact5 > wtor=weights(13)*fact > wtor_d=weights(14)*fact2 > wsccor=weights(21)*fact > 237c901 < subroutine enerprint(energia,fact) --- > subroutine enerprint(energia) 240d903 < include 'DIMENSIONS.ZSCOPT' 244c907,908 < double precision energia(0:max_ene),fact(6) --- > include 'COMMON.MD' > double precision energia(0:n_ene) 246c910,911 < evdw=energia(1)+fact(6)*energia(21) --- > evdw=energia(1) > evdw2=energia(2) 248c913 < evdw2=energia(2)+energia(17) --- > evdw2=energia(2)+energia(18) 268,270c933,936 < esccor=energia(19) < edihcnstr=energia(20) < estr=energia(18) --- > edihcnstr=energia(19) > estr=energia(17) > Uconst=energia(20) > esccor=energia(21) 272,279c938,945 < write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, < & wvdwpp, < & estr,wbond,ebe,wang,escloc,wscloc,etors,wtor*fact(1), < & etors_d,wtor_d*fact(2),ehpb,wstrain, < & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), < & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), < & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), < & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot --- > write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, > & estr,wbond,ebe,wang, > & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, > & ecorr,wcorr, > & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, > & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, > & edihcnstr,ebr*nss, > & Uconst,etot 283c949 < & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p elec)'/ --- > & 'EES= ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/ 301c967,968 < & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ --- > & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ > & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 304,310c971,977 < write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond, < & ebe,wang,escloc,wscloc,etors,wtor*fact(1),etors_d,wtor_d*fact2, < & ehpb,wstrain,ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4), < & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), < & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), < & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, < & edihcnstr,ebr*nss,etot --- > write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, > & estr,wbond,ebe,wang, > & escloc,wscloc,etors,wtor,etors_d,wtor_d,ehpb,wstrain, > & ecorr,wcorr, > & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, > & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, > & ebr*nss,Uconst,etot 331c998,999 < & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ --- > & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ > & 'UCONST=',1pE16.6,' (Constraint energy)'/ 337c1005 < subroutine elj(evdw,evdw_t) --- > subroutine elj(evdw) 344,345d1011 < include 'DIMENSIONS.ZSCOPT' < include "DIMENSIONS.COMPAR" 354d1019 < include 'COMMON.ENEPS' 360,367c1025 < integer icant < external icant < cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon < do i=1,210 < do j=1,2 < eneps_temp(j,i)=0.0d0 < enddo < enddo --- > c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon 369d1026 < evdw_t=0.0d0 400,402d1056 < ij=icant(itypi,itypj) < eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) < eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij 409,414c1063 < if (bb(itypi,itypj).gt.0.0d0) then < evdw=evdw+evdwij < else < evdw_t=evdw_t+evdwij < endif < if (calc_grad) then --- > evdw=evdw+evdwij 424a1074,1075 > gvdwc(k,i)=gvdwc(k,i)-gg(k) > gvdwc(k,j)=gvdwc(k,j)+gg(k) 426,431c1077,1081 < do k=i,j-1 < do l=1,3 < gvdwc(l,k)=gvdwc(l,k)+gg(l) < enddo < enddo < endif --- > cgrad do k=i,j-1 > cgrad do l=1,3 > cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) > cgrad enddo > cgrad enddo 493d1142 < if (calc_grad) then 500d1148 < endif 513c1161 < subroutine eljk(evdw,evdw_t) --- > subroutine eljk(evdw) 520,521d1167 < include 'DIMENSIONS.ZSCOPT' < include "DIMENSIONS.COMPAR" 528d1173 < include 'COMMON.ENEPS' 533,534d1177 < integer icant < external icant 536,540d1178 < do i=1,210 < do j=1,2 < eneps_temp(j,i)=0.0d0 < enddo < enddo 542d1179 < evdw_t=0.0d0 570,573d1206 < ij=icant(itypi,itypj) < eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm) < & /dabs(eps(itypi,itypj)) < eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps(itypi,itypj) 581,586c1214 < if (bb(itypi,itypj).gt.0.0d0) then < evdw=evdw+evdwij < else < evdw_t=evdw_t+evdwij < endif < if (calc_grad) then --- > evdw=evdw+evdwij 596a1225,1226 > gvdwc(k,i)=gvdwc(k,i)-gg(k) > gvdwc(k,j)=gvdwc(k,j)+gg(k) 598,603c1228,1232 < do k=i,j-1 < do l=1,3 < gvdwc(l,k)=gvdwc(l,k)+gg(l) < enddo < enddo < endif --- > cgrad do k=i,j-1 > cgrad do l=1,3 > cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) > cgrad enddo > cgrad enddo 607d1235 < if (calc_grad) then 614d1241 < endif 618c1245 < subroutine ebp(evdw,evdw_t) --- > subroutine ebp(evdw) 625,626d1251 < include 'DIMENSIONS.ZSCOPT' < include "DIMENSIONS.COMPAR" 634d1258 < include 'COMMON.ENEPS' 640,646d1263 < integer icant < external icant < do i=1,210 < do j=1,2 < eneps_temp(j,i)=0.0d0 < enddo < enddo 648d1264 < evdw_t=0.0d0 649a1266 > evdw=0.0D0 665a1283 > c dsci_inv=dsc_inv(itypi) 674a1293 > c dscj_inv=dsc_inv(itypj) 719,729c1338 < ij=icant(itypi,itypj) < aux=eps1*eps2rt**2*eps3rt**2 < eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux < & /dabs(eps(itypi,itypj)) < eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj) < if (bb(itypi,itypj).gt.0.0d0) then < evdw=evdw+evdwij < else < evdw_t=evdw_t+evdwij < endif < if (calc_grad) then --- > evdw=evdw+evdwij 752d1360 < endif 760c1368 < subroutine egb(evdw,evdw_t) --- > subroutine egb(evdw) 767,768d1374 < include 'DIMENSIONS.ZSCOPT' < include "DIMENSIONS.COMPAR" 776d1381 < include 'COMMON.ENEPS' 778a1384 > include 'COMMON.CONTROL' 780,787c1386,1387 < common /srutu/icall < integer icant < external icant < do i=1,210 < do j=1,2 < eneps_temp(j,i)=0.0d0 < enddo < enddo --- > evdw=0.0D0 > ccccc energy_dec=.false. 790d1389 < evdw_t=0.0d0 792c1391 < c if (icall.gt.0) lprn=.true. --- > c if (icall.eq.0) lprn=.false. 803a1403 > c dsci_inv=dsc_inv(itypi) 804a1405,1406 > c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres) > c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi 812a1415 > c dscj_inv=dsc_inv(itypj) 813a1417,1419 > c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, > c & 1.0d0/vbld(j+nres) > c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 840c1446,1448 < c write (iout,*) i,j,xj,yj,zj --- > c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi > c write (iout,*) "j",j," dc_norm", > c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j) 848a1457,1458 > c for diagnostics; uncomment > c rij_shift=1.2*sig0ij 851a1462,1464 > cd write (iout,'(2(a3,i3,2x),17(0pf7.3))') > cd & restyp(itypi),i,restyp(itypj),j, > cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 862a1476,1477 > c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt, > c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 864,876c1479 < if (bb(itypi,itypj).gt.0) then < evdw=evdw+evdwij < else < evdw_t=evdw_t+evdwij < endif < ij=icant(itypi,itypj) < aux=eps1*eps2rt**2*eps3rt**2 < eneps_temp(1,ij)=eneps_temp(1,ij)+aux*e1 < & /dabs(eps(itypi,itypj)) < eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj) < c write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj, < c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)), < c & aux*e2/eps(itypi,itypj) --- > evdw=evdw+evdwij 887c1490,1493 < if (calc_grad) then --- > > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'evdw',i,j,evdwij > 892a1499 > c fac=0.0d0 899d1505 < endif 902a1509,1510 > c write (iout,*) "Number of loop steps in EGB:",ind > cccc energy_dec=.false. 906c1514 < subroutine egbv(evdw,evdw_t) --- > subroutine egbv(evdw) 913,914d1520 < include 'DIMENSIONS.ZSCOPT' < include "DIMENSIONS.COMPAR" 922d1527 < include 'COMMON.ENEPS' 927,933d1531 < integer icant < external icant < do i=1,210 < do j=1,2 < eneps_temp(j,i)=0.0d0 < enddo < enddo 935d1532 < evdw_t=0.0d0 939c1536 < c if (icall.gt.0) lprn=.true. --- > c if (icall.eq.0) lprn=.true. 950a1548 > c dsci_inv=dsc_inv(itypi) 959a1558 > c dscj_inv=dsc_inv(itypj) 1013,1016c1612,1622 < if (bb(itypi,itypj).gt.0.0d0) then < evdw=evdw+evdwij+e_augm < else < evdw_t=evdw_t+evdwij+e_augm --- > evdw=evdw+evdwij+e_augm > if (lprn) then > sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) > epsi=bb(itypi,itypj)**2/aa(itypi,itypj) > write (iout,'(2(a3,i3,2x),17(0pf7.3))') > & restyp(itypi),i,restyp(itypj),j, > & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0), > & chi1,chi2,chip1,chip2, > & eps1,eps2rt**2,eps3rt**2, > & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, > & evdwij+e_augm 1018,1036d1623 < ij=icant(itypi,itypj) < aux=eps1*eps2rt**2*eps3rt**2 < eneps_temp(1,ij)=eneps_temp(1,ij)+aux*(e1+e_augm) < & /dabs(eps(itypi,itypj)) < eneps_temp(2,ij)=eneps_temp(2,ij)+aux*e2/eps(itypi,itypj) < c eneps_temp(ij)=eneps_temp(ij) < c & +(evdwij+e_augm)/eps(itypi,itypj) < c if (lprn) then < c sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) < c epsi=bb(itypi,itypj)**2/aa(itypi,itypj) < c write (iout,'(2(a3,i3,2x),17(0pf7.3))') < c & restyp(itypi),i,restyp(itypj),j, < c & epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0), < c & chi1,chi2,chip1,chip2, < c & eps1,eps2rt**2,eps3rt**2, < c & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, < c & evdwij+e_augm < c endif < if (calc_grad) then 1048d1634 < endif 1052d1637 < return 1059a1645 > include 'COMMON.IOUNITS' 1072a1659,1663 > c diagnostics only > c faceps1_inv=om12 > c eps1=om12 > c eps1_om12=1.0d0 > c write (iout,*) "om12",om12," eps1",eps1 1082a1674,1681 > c diagnostics only > c sigsq=1.0d0 > c sigsq_om1=0.0d0 > c sigsq_om2=0.0d0 > c sigsq_om12=0.0d0 > c write (iout,*) "chiom1",chiom1," chiom2",chiom2," chiom12",chiom12 > c write (iout,*) "faceps1",faceps1," faceps1_inv",faceps1_inv, > c & " eps1",eps1 1089a1689,1690 > c write (iout,*) "chipom1",chipom1," chipom2",chipom2, > c & " chipom12",chipom12," facp",facp," facp_inv",facp_inv 1098a1700,1702 > c write (iout,*) "eps2rt",eps2rt," eps3rt",eps3rt > c write (iout,*) "eps2rt_om1",eps2rt_om1," eps2rt_om2",eps2rt_om2, > c & " eps2rt_om12",eps2rt_om12 1107d1710 < include 'DIMENSIONS.ZSCOPT' 1110a1714 > include 'COMMON.IOUNITS' 1115a1720,1728 > c diagnostics only > c eom1=0.0d0 > c eom2=0.0d0 > c eom12=evdwij*eps1_om12 > c end diagnostics > c write (iout,*) "eps2der",eps2der," eps3der",eps3der, > c & " sigder",sigder > c write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12 > c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 1122a1736 > c write (iout,*) "gg",(gg(k),k=1,3) 1129a1744,1747 > c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) > c & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv > c write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) > c & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv 1134,1137c1752,1759 < do k=i,j-1 < do l=1,3 < gvdwc(l,k)=gvdwc(l,k)+gg(l) < enddo --- > cgrad do k=i,j-1 > cgrad do l=1,3 > cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) > cgrad enddo > cgrad enddo > do l=1,3 > gvdwc(l,i)=gvdwc(l,i)-gg(l) > gvdwc(l,j)=gvdwc(l,j)+gg(l) 1141,1142c1763,1768 < c------------------------------------------------------------------------------ < subroutine vec_and_deriv --- > C----------------------------------------------------------------------- > subroutine e_softsphere(evdw) > C > C This subroutine calculates the interaction energy of nonbonded side chains > C assuming the LJ potential of interaction. > C 1145,1146c1771 < include 'DIMENSIONS.ZSCOPT' < include 'COMMON.IOUNITS' --- > parameter (accur=1.0d-10) 1151d1775 < include 'COMMON.VECTORS' 1154,1247c1778,1815 < dimension uyder(3,3,2),uzder(3,3,2),vbld_inv_temp(2) < C Compute the local reference systems. For reference system (i), the < C X-axis points from CA(i) to CA(i+1), the Y axis is in the < C CA(i)-CA(i+1)-CA(i+2) plane, and the Z axis is perpendicular to this plane. < do i=1,nres-1 < c if (i.eq.nres-1 .or. itel(i+1).eq.0) then < if (i.eq.nres-1) then < C Case of the last full residue < C Compute the Z-axis < call vecpr(dc_norm(1,i),dc_norm(1,i-1),uz(1,i)) < costh=dcos(pi-theta(nres)) < fac=1.0d0/dsqrt(1.0d0-costh*costh) < do k=1,3 < uz(k,i)=fac*uz(k,i) < enddo < if (calc_grad) then < C Compute the derivatives of uz < uzder(1,1,1)= 0.0d0 < uzder(2,1,1)=-dc_norm(3,i-1) < uzder(3,1,1)= dc_norm(2,i-1) < uzder(1,2,1)= dc_norm(3,i-1) < uzder(2,2,1)= 0.0d0 < uzder(3,2,1)=-dc_norm(1,i-1) < uzder(1,3,1)=-dc_norm(2,i-1) < uzder(2,3,1)= dc_norm(1,i-1) < uzder(3,3,1)= 0.0d0 < uzder(1,1,2)= 0.0d0 < uzder(2,1,2)= dc_norm(3,i) < uzder(3,1,2)=-dc_norm(2,i) < uzder(1,2,2)=-dc_norm(3,i) < uzder(2,2,2)= 0.0d0 < uzder(3,2,2)= dc_norm(1,i) < uzder(1,3,2)= dc_norm(2,i) < uzder(2,3,2)=-dc_norm(1,i) < uzder(3,3,2)= 0.0d0 < endif < C Compute the Y-axis < facy=fac < do k=1,3 < uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i)) < enddo < if (calc_grad) then < C Compute the derivatives of uy < do j=1,3 < do k=1,3 < uyder(k,j,1)=2*dc_norm(k,i-1)*dc_norm(j,i) < & -dc_norm(k,i)*dc_norm(j,i-1) < uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i) < enddo < uyder(j,j,1)=uyder(j,j,1)-costh < uyder(j,j,2)=1.0d0+uyder(j,j,2) < enddo < do j=1,2 < do k=1,3 < do l=1,3 < uygrad(l,k,j,i)=uyder(l,k,j) < uzgrad(l,k,j,i)=uzder(l,k,j) < enddo < enddo < enddo < call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i)) < call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i)) < call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i)) < call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i)) < endif < else < C Other residues < C Compute the Z-axis < call vecpr(dc_norm(1,i),dc_norm(1,i+1),uz(1,i)) < costh=dcos(pi-theta(i+2)) < fac=1.0d0/dsqrt(1.0d0-costh*costh) < do k=1,3 < uz(k,i)=fac*uz(k,i) < enddo < if (calc_grad) then < C Compute the derivatives of uz < uzder(1,1,1)= 0.0d0 < uzder(2,1,1)=-dc_norm(3,i+1) < uzder(3,1,1)= dc_norm(2,i+1) < uzder(1,2,1)= dc_norm(3,i+1) < uzder(2,2,1)= 0.0d0 < uzder(3,2,1)=-dc_norm(1,i+1) < uzder(1,3,1)=-dc_norm(2,i+1) < uzder(2,3,1)= dc_norm(1,i+1) < uzder(3,3,1)= 0.0d0 < uzder(1,1,2)= 0.0d0 < uzder(2,1,2)= dc_norm(3,i) < uzder(3,1,2)=-dc_norm(2,i) < uzder(1,2,2)=-dc_norm(3,i) < uzder(2,2,2)= 0.0d0 < uzder(3,2,2)= dc_norm(1,i) < uzder(1,3,2)= dc_norm(2,i) < uzder(2,3,2)=-dc_norm(1,i) < uzder(3,3,2)= 0.0d0 --- > include 'COMMON.TORSION' > include 'COMMON.SBRIDGE' > include 'COMMON.NAMES' > include 'COMMON.IOUNITS' > include 'COMMON.CONTACTS' > dimension gg(3) > cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct > evdw=0.0D0 > do i=iatsc_s,iatsc_e > itypi=itype(i) > if (itypi.eq.21) cycle > itypi1=itype(i+1) > xi=c(1,nres+i) > yi=c(2,nres+i) > zi=c(3,nres+i) > C > C Calculate SC interaction energy. > C > do iint=1,nint_gr(i) > cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), > cd & 'iend=',iend(i,iint) > do j=istart(i,iint),iend(i,iint) > itypj=itype(j) > if (itypj.eq.21) cycle > xj=c(1,nres+j)-xi > yj=c(2,nres+j)-yi > zj=c(3,nres+j)-zi > rij=xj*xj+yj*yj+zj*zj > c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj > r0ij=r0(itypi,itypj) > r0ijsq=r0ij*r0ij > c print *,i,j,r0ij,dsqrt(rij) > if (rij.lt.r0ijsq) then > evdwij=0.25d0*(rij-r0ijsq)**2 > fac=rij-r0ijsq > else > evdwij=0.0d0 > fac=0.0d0 1249,1250c1817,1823 < C Compute the Y-axis < facy=fac --- > evdw=evdw+evdwij > C > C Calculate the components of the gradient in DC and X > C > gg(1)=xj*fac > gg(2)=yj*fac > gg(3)=zj*fac 1252,1263c1825,1828 < uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i)) < enddo < if (calc_grad) then < C Compute the derivatives of uy < do j=1,3 < do k=1,3 < uyder(k,j,1)=2*dc_norm(k,i+1)*dc_norm(j,i) < & -dc_norm(k,i)*dc_norm(j,i+1) < uyder(k,j,2)=-dc_norm(j,i)*dc_norm(k,i) < enddo < uyder(j,j,1)=uyder(j,j,1)-costh < uyder(j,j,2)=1.0d0+uyder(j,j,2) --- > gvdwx(k,i)=gvdwx(k,i)-gg(k) > gvdwx(k,j)=gvdwx(k,j)+gg(k) > gvdwc(k,i)=gvdwc(k,i)-gg(k) > gvdwc(k,j)=gvdwc(k,j)+gg(k) 1265,1277c1830,1898 < do j=1,2 < do k=1,3 < do l=1,3 < uygrad(l,k,j,i)=uyder(l,k,j) < uzgrad(l,k,j,i)=uzder(l,k,j) < enddo < enddo < enddo < call unormderiv(uy(1,i),uyder(1,1,1),facy,uygrad(1,1,1,i)) < call unormderiv(uy(1,i),uyder(1,1,2),facy,uygrad(1,1,2,i)) < call unormderiv(uz(1,i),uzder(1,1,1),fac,uzgrad(1,1,1,i)) < call unormderiv(uz(1,i),uzder(1,1,2),fac,uzgrad(1,1,2,i)) < endif --- > cgrad do k=i,j-1 > cgrad do l=1,3 > cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) > cgrad enddo > cgrad enddo > enddo ! j > enddo ! iint > enddo ! i > return > end > C-------------------------------------------------------------------------- > subroutine eelec_soft_sphere(ees,evdw1,eel_loc,eello_turn3, > & eello_turn4) > C > C Soft-sphere potential of p-p interaction > C > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'COMMON.CONTROL' > include 'COMMON.IOUNITS' > include 'COMMON.GEO' > include 'COMMON.VAR' > include 'COMMON.LOCAL' > include 'COMMON.CHAIN' > include 'COMMON.DERIV' > include 'COMMON.INTERACT' > include 'COMMON.CONTACTS' > include 'COMMON.TORSION' > include 'COMMON.VECTORS' > include 'COMMON.FFIELD' > dimension ggg(3) > cd write(iout,*) 'In EELEC_soft_sphere' > ees=0.0D0 > evdw1=0.0D0 > eel_loc=0.0d0 > eello_turn3=0.0d0 > eello_turn4=0.0d0 > ind=0 > do i=iatel_s,iatel_e > if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle > dxi=dc(1,i) > dyi=dc(2,i) > dzi=dc(3,i) > xmedi=c(1,i)+0.5d0*dxi > ymedi=c(2,i)+0.5d0*dyi > zmedi=c(3,i)+0.5d0*dzi > num_conti=0 > c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) > do j=ielstart(i),ielend(i) > if (itype(j).eq.21 .or. itype(j+1).eq.21) cycle > ind=ind+1 > iteli=itel(i) > itelj=itel(j) > if (j.eq.i+2 .and. itelj.eq.2) iteli=2 > r0ij=rpp(iteli,itelj) > r0ijsq=r0ij*r0ij > dxj=dc(1,j) > dyj=dc(2,j) > dzj=dc(3,j) > xj=c(1,j)+0.5D0*dxj-xmedi > yj=c(2,j)+0.5D0*dyj-ymedi > zj=c(3,j)+0.5D0*dzj-zmedi > rij=xj*xj+yj*yj+zj*zj > if (rij.lt.r0ijsq) then > evdw1ij=0.25d0*(rij-r0ijsq)**2 > fac=rij-r0ijsq > else > evdw1ij=0.0d0 > fac=0.0d0 1279,1288c1900,1906 < enddo < if (calc_grad) then < do i=1,nres-1 < vbld_inv_temp(1)=vbld_inv(i+1) < if (i.lt.nres-1) then < vbld_inv_temp(2)=vbld_inv(i+2) < else < vbld_inv_temp(2)=vbld_inv(i) < endif < do j=1,2 --- > evdw1=evdw1+evdw1ij > C > C Calculate contributions to the Cartesian gradient. > C > ggg(1)=fac*xj > ggg(2)=fac*yj > ggg(3)=fac*zj 1290,1293c1908,1909 < do l=1,3 < uygrad(l,k,j,i)=vbld_inv_temp(j)*uygrad(l,k,j,i) < uzgrad(l,k,j,i)=vbld_inv_temp(j)*uzgrad(l,k,j,i) < enddo --- > gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) > gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) 1295,1297c1911,1930 < enddo < enddo < endif --- > * > * Loop over residues i+1 thru j-1. > * > cgrad do k=i+1,j-1 > cgrad do l=1,3 > cgrad gelc(l,k)=gelc(l,k)+ggg(l) > cgrad enddo > cgrad enddo > enddo ! j > enddo ! i > cgrad do i=nnt,nct-1 > cgrad do k=1,3 > cgrad gelc(k,i)=gelc(k,i)+0.5d0*gelc(k,i) > cgrad enddo > cgrad do j=i+1,nct-1 > cgrad do k=1,3 > cgrad gelc(k,i)=gelc(k,i)+gelc(k,j) > cgrad enddo > cgrad enddo > cgrad enddo 1300,1301c1933,1934 < C----------------------------------------------------------------------------- < subroutine vec_and_deriv_test --- > c------------------------------------------------------------------------------ > subroutine vec_and_deriv 1304c1937,1939 < include 'DIMENSIONS.ZSCOPT' --- > #ifdef MPI > include 'mpif.h' > #endif 1311c1946,1948 < dimension uyder(3,3,2),uzder(3,3,2) --- > include 'COMMON.SETUP' > include 'COMMON.TIME1' > dimension uyder(3,3,2),uzder(3,3,2),vbld_inv_temp(2) 1314a1952,1954 > #ifdef PARVEC > do i=ivec_start,ivec_end > #else 1315a1956 > #endif 1322,1324d1962 < c write (iout,*) 'fac',fac, < c & 1.0d0/dsqrt(scalar(uz(1,i),uz(1,i))) < fac=1.0d0/dsqrt(scalar(uz(1,i),uz(1,i))) 1348,1350d1985 < do k=1,3 < uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i)) < enddo 1352,1365d1986 < facy=1.0d0/dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))* < & (scalar(dc_norm(1,i-1),dc_norm(1,i-1))**2- < & scalar(dc_norm(1,i),dc_norm(1,i-1))**2)) < do k=1,3 < c uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i)) < uy(k,i)= < c & facy*( < & dc_norm(k,i-1)*scalar(dc_norm(1,i),dc_norm(1,i)) < & -scalar(dc_norm(1,i),dc_norm(1,i-1))*dc_norm(k,i) < c & ) < enddo < c write (iout,*) 'facy',facy, < c & 1.0d0/dsqrt(scalar(uy(1,i),uy(1,i))) < facy=1.0d0/dsqrt(scalar(uy(1,i),uy(1,i))) 1367c1988 < uy(k,i)=facy*uy(k,i) --- > uy(k,i)=fac*(dc_norm(k,i-1)-costh*dc_norm(k,i)) 1376,1381c1997,1998 < c uyder(j,j,1)=uyder(j,j,1)-costh < c uyder(j,j,2)=1.0d0+uyder(j,j,2) < uyder(j,j,1)=uyder(j,j,1) < & -scalar(dc_norm(1,i),dc_norm(1,i-1)) < uyder(j,j,2)=scalar(dc_norm(1,i),dc_norm(1,i)) < & +uyder(j,j,2) --- > uyder(j,j,1)=uyder(j,j,1)-costh > uyder(j,j,2)=1.0d0+uyder(j,j,2) 1401d2017 < fac=1.0d0/dsqrt(scalar(uz(1,i),uz(1,i))) 1426,1439d2041 < facy=1.0d0/dsqrt(scalar(dc_norm(1,i),dc_norm(1,i))* < & (scalar(dc_norm(1,i+1),dc_norm(1,i+1))**2- < & scalar(dc_norm(1,i),dc_norm(1,i+1))**2)) < do k=1,3 < c uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i)) < uy(k,i)= < c & facy*( < & dc_norm(k,i+1)*scalar(dc_norm(1,i),dc_norm(1,i)) < & -scalar(dc_norm(1,i),dc_norm(1,i+1))*dc_norm(k,i) < c & ) < enddo < c write (iout,*) 'facy',facy, < c & 1.0d0/dsqrt(scalar(uy(1,i),uy(1,i))) < facy=1.0d0/dsqrt(scalar(uy(1,i),uy(1,i))) 1441c2043 < uy(k,i)=facy*uy(k,i) --- > uy(k,i)=facy*(dc_norm(k,i+1)-costh*dc_norm(k,i)) 1450,1455c2052,2053 < c uyder(j,j,1)=uyder(j,j,1)-costh < c uyder(j,j,2)=1.0d0+uyder(j,j,2) < uyder(j,j,1)=uyder(j,j,1) < & -scalar(dc_norm(1,i),dc_norm(1,i+1)) < uyder(j,j,2)=scalar(dc_norm(1,i),dc_norm(1,i)) < & +uyder(j,j,2) --- > uyder(j,j,1)=uyder(j,j,1)-costh > uyder(j,j,2)=1.0d0+uyder(j,j,2) 1471a2070,2075 > vbld_inv_temp(1)=vbld_inv(i+1) > if (i.lt.nres-1) then > vbld_inv_temp(2)=vbld_inv(i+2) > else > vbld_inv_temp(2)=vbld_inv(i) > endif 1475,1476c2079,2080 < uygrad(l,k,j,i)=vblinv*uygrad(l,k,j,i) < uzgrad(l,k,j,i)=vblinv*uzgrad(l,k,j,i) --- > uygrad(l,k,j,i)=vbld_inv_temp(j)*uygrad(l,k,j,i) > uzgrad(l,k,j,i)=vbld_inv_temp(j)*uzgrad(l,k,j,i) 1480a2085,2112 > #if defined(PARVEC) && defined(MPI) > if (nfgtasks1.gt.1) then > time00=MPI_Wtime() > c print *,"Processor",fg_rank1,kolor1," ivec_start",ivec_start, > c & " ivec_displ",(ivec_displ(i),i=0,nfgtasks1-1), > c & " ivec_count",(ivec_count(i),i=0,nfgtasks1-1) > call MPI_Allgatherv(uy(1,ivec_start),ivec_count(fg_rank1), > & MPI_UYZ,uy(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ, > & FG_COMM1,IERR) > call MPI_Allgatherv(uz(1,ivec_start),ivec_count(fg_rank1), > & MPI_UYZ,uz(1,1),ivec_count(0),ivec_displ(0),MPI_UYZ, > & FG_COMM1,IERR) > call MPI_Allgatherv(uygrad(1,1,1,ivec_start), > & ivec_count(fg_rank1),MPI_UYZGRAD,uygrad(1,1,1,1),ivec_count(0), > & ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR) > call MPI_Allgatherv(uzgrad(1,1,1,ivec_start), > & ivec_count(fg_rank1),MPI_UYZGRAD,uzgrad(1,1,1,1),ivec_count(0), > & ivec_displ(0),MPI_UYZGRAD,FG_COMM1,IERR) > time_gather=time_gather+MPI_Wtime()-time00 > endif > c if (fg_rank.eq.0) then > c write (iout,*) "Arrays UY and UZ" > c do i=1,nres-1 > c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(uy(k,i),k=1,3), > c & (uz(k,i),k=1,3) > c enddo > c endif > #endif 1487d2118 < include 'DIMENSIONS.ZSCOPT' 1572c2203,2208 < include 'DIMENSIONS.ZSCOPT' --- > #ifdef MPI > include "mpif.h" > include "COMMON.SETUP" > integer IERR > integer status(MPI_STATUS_SIZE) > #endif 1588a2225,2227 > #ifdef PARMAT > do i=ivec_start+2,ivec_end+2 > #else 1589a2229 > #endif 1655a2296 > c if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then 1657,1661c2298 < if (itype(i-2).le.ntyp) then < iti = itortyp(itype(i-2)) < else < iti=ntortyp+1 < endif --- > iti = itortyp(itype(i-2)) 1664a2302 > c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then 1666,1670c2304 < if (itype(i-1).le.ntyp) then < iti1 = itortyp(itype(i-1)) < else < iti1=ntortyp+1 < endif --- > iti1 = itortyp(itype(i-1)) 1678,1679c2312,2313 < c print *,"itilde1 i iti iti1",i,iti,iti1 < if (i .gt. iatel_s+2) then --- > c if (i .gt. iatel_s+2) then > if (i .gt. nnt+2) then 1681a2316,2317 > if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) > & then 1686a2323 > endif 1700d2336 < c print *,"itilde2 i iti iti1",i,iti,iti1 1703,1708d2338 < call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2)) < call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2)) < call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2)) < call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2)) < call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2)) < c print *,"itilde3 i iti iti1",i,iti,iti1 1711a2342 > c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then 1713,1717c2344 < if (itype(i-1).le.ntyp) then < iti1 = itortyp(itype(i-1)) < else < iti1=ntortyp+1 < endif --- > iti1 = itortyp(itype(i-1)) 1723a2351,2360 > cd write (iout,*) 'mu ',mu(:,i-2) > cd write (iout,*) 'mu1',mu1(:,i-2) > cd write (iout,*) 'mu2',mu2(:,i-2) > if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) > & then > call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2)) > call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2)) > call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2)) > call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2)) > call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2)) 1734,1735c2371 < cd write (iout,*) 'i',i,' mu ',(mu(k,i-2),k=1,2), < cd & ' mu1',(b1(k,i-2),k=1,2),' mu2',(Ub2(k,i-2),k=1,2) --- > endif 1738a2375,2377 > if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) > &then > c do i=max0(ivec_start,2),ivec_end 1748a2388,2636 > endif > #if defined(MPI) && defined(PARMAT) > #ifdef DEBUG > c if (fg_rank.eq.0) then > write (iout,*) "Arrays UG and UGDER before GATHER" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & ((ug(l,k,i),l=1,2),k=1,2), > & ((ugder(l,k,i),l=1,2),k=1,2) > enddo > write (iout,*) "Arrays UG2 and UG2DER" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & ((ug2(l,k,i),l=1,2),k=1,2), > & ((ug2der(l,k,i),l=1,2),k=1,2) > enddo > write (iout,*) "Arrays OBROT OBROT2 OBROTDER and OBROT2DER" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & (obrot(k,i),k=1,2),(obrot2(k,i),k=1,2), > & (obrot_der(k,i),k=1,2),(obrot2_der(k,i),k=1,2) > enddo > write (iout,*) "Arrays COSTAB SINTAB COSTAB2 and SINTAB2" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & costab(i),sintab(i),costab2(i),sintab2(i) > enddo > write (iout,*) "Array MUDER" > do i=1,nres-1 > write (iout,'(i5,2f10.5)') i,muder(1,i),muder(2,i) > enddo > c endif > #endif > if (nfgtasks.gt.1) then > time00=MPI_Wtime() > c write(iout,*)"Processor",fg_rank,kolor," ivec_start",ivec_start, > c & " ivec_displ",(ivec_displ(i),i=0,nfgtasks-1), > c & " ivec_count",(ivec_count(i),i=0,nfgtasks-1) > #ifdef MATGATHER > call MPI_Allgatherv(Ub2(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Ub2(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Ub2der(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Ub2der(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(mu(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,mu(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(muder(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,muder(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Eug(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Eug(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Eugder(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Eugder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(costab(ivec_start),ivec_count(fg_rank1), > & MPI_DOUBLE_PRECISION,costab(1),ivec_count(0),ivec_displ(0), > & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) > call MPI_Allgatherv(sintab(ivec_start),ivec_count(fg_rank1), > & MPI_DOUBLE_PRECISION,sintab(1),ivec_count(0),ivec_displ(0), > & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) > call MPI_Allgatherv(costab2(ivec_start),ivec_count(fg_rank1), > & MPI_DOUBLE_PRECISION,costab2(1),ivec_count(0),ivec_displ(0), > & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) > call MPI_Allgatherv(sintab2(ivec_start),ivec_count(fg_rank1), > & MPI_DOUBLE_PRECISION,sintab2(1),ivec_count(0),ivec_displ(0), > & MPI_DOUBLE_PRECISION,FG_COMM1,IERR) > if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) > & then > call MPI_Allgatherv(Ctobr(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Ctobr(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Ctobrder(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Ctobrder(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Dtobr2(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Dtobr2(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Dtobr2der(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Dtobr2der(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Ug2Db1t(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,Ug2Db1t(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Ug2Db1tder(1,ivec_start), > & ivec_count(fg_rank1), > & MPI_MU,Ug2Db1tder(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(CUgb2(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,CUgb2(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(CUgb2der(1,ivec_start),ivec_count(fg_rank1), > & MPI_MU,CUgb2der(1,1),ivec_count(0),ivec_displ(0),MPI_MU, > & FG_COMM1,IERR) > call MPI_Allgatherv(Cug(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Cug(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Cugder(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Cugder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Dug(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Dug(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Dugder(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Dugder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Dtug2(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,Dtug2(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Dtug2der(1,1,ivec_start), > & ivec_count(fg_rank1), > & MPI_MAT1,Dtug2der(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(EugC(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,EugC(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(EugCder(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,EugCder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(EugD(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,EugD(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(EugDder(1,1,ivec_start),ivec_count(fg_rank1), > & MPI_MAT1,EugDder(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(DtUg2EUg(1,1,ivec_start), > & ivec_count(fg_rank1), > & MPI_MAT1,DtUg2EUg(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(Ug2DtEUg(1,1,ivec_start), > & ivec_count(fg_rank1), > & MPI_MAT1,Ug2DtEUg(1,1,1),ivec_count(0),ivec_displ(0),MPI_MAT1, > & FG_COMM1,IERR) > call MPI_Allgatherv(DtUg2EUgder(1,1,1,ivec_start), > & ivec_count(fg_rank1), > & MPI_MAT2,DtUg2EUgder(1,1,1,1),ivec_count(0),ivec_displ(0), > & MPI_MAT2,FG_COMM1,IERR) > call MPI_Allgatherv(Ug2DtEUgder(1,1,1,ivec_start), > & ivec_count(fg_rank1), > & MPI_MAT2,Ug2DtEUgder(1,1,1,1),ivec_count(0),ivec_displ(0), > & MPI_MAT2,FG_COMM1,IERR) > endif > #else > c Passes matrix info through the ring > isend=fg_rank1 > irecv=fg_rank1-1 > if (irecv.lt.0) irecv=nfgtasks1-1 > iprev=irecv > inext=fg_rank1+1 > if (inext.ge.nfgtasks1) inext=0 > do i=1,nfgtasks1-1 > c write (iout,*) "isend",isend," irecv",irecv > c call flush(iout) > lensend=lentyp(isend) > lenrecv=lentyp(irecv) > c write (iout,*) "lensend",lensend," lenrecv",lenrecv > c call MPI_SENDRECV(ug(1,1,ivec_displ(isend)+1),1, > c & MPI_ROTAT1(lensend),inext,2200+isend, > c & ug(1,1,ivec_displ(irecv)+1),1,MPI_ROTAT1(lenrecv), > c & iprev,2200+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather ROTAT1" > c call flush(iout) > c call MPI_SENDRECV(obrot(1,ivec_displ(isend)+1),1, > c & MPI_ROTAT2(lensend),inext,3300+isend, > c & obrot(1,ivec_displ(irecv)+1),1,MPI_ROTAT2(lenrecv), > c & iprev,3300+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather ROTAT2" > c call flush(iout) > call MPI_SENDRECV(costab(ivec_displ(isend)+1),1, > & MPI_ROTAT_OLD(lensend),inext,4400+isend, > & costab(ivec_displ(irecv)+1),1,MPI_ROTAT_OLD(lenrecv), > & iprev,4400+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather ROTAT_OLD" > c call flush(iout) > call MPI_SENDRECV(mu(1,ivec_displ(isend)+1),1, > & MPI_PRECOMP11(lensend),inext,5500+isend, > & mu(1,ivec_displ(irecv)+1),1,MPI_PRECOMP11(lenrecv), > & iprev,5500+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather PRECOMP11" > c call flush(iout) > call MPI_SENDRECV(Eug(1,1,ivec_displ(isend)+1),1, > & MPI_PRECOMP12(lensend),inext,6600+isend, > & Eug(1,1,ivec_displ(irecv)+1),1,MPI_PRECOMP12(lenrecv), > & iprev,6600+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather PRECOMP12" > c call flush(iout) > if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) > & then > call MPI_SENDRECV(ug2db1t(1,ivec_displ(isend)+1),1, > & MPI_ROTAT2(lensend),inext,7700+isend, > & ug2db1t(1,ivec_displ(irecv)+1),1,MPI_ROTAT2(lenrecv), > & iprev,7700+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather PRECOMP21" > c call flush(iout) > call MPI_SENDRECV(EUgC(1,1,ivec_displ(isend)+1),1, > & MPI_PRECOMP22(lensend),inext,8800+isend, > & EUgC(1,1,ivec_displ(irecv)+1),1,MPI_PRECOMP22(lenrecv), > & iprev,8800+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather PRECOMP22" > c call flush(iout) > call MPI_SENDRECV(Ug2DtEUgder(1,1,1,ivec_displ(isend)+1),1, > & MPI_PRECOMP23(lensend),inext,9900+isend, > & Ug2DtEUgder(1,1,1,ivec_displ(irecv)+1),1, > & MPI_PRECOMP23(lenrecv), > & iprev,9900+irecv,FG_COMM,status,IERR) > c write (iout,*) "Gather PRECOMP23" > c call flush(iout) > endif > isend=irecv > irecv=irecv-1 > if (irecv.lt.0) irecv=nfgtasks1-1 > enddo > #endif > time_gather=time_gather+MPI_Wtime()-time00 > endif > #ifdef DEBUG > c if (fg_rank.eq.0) then > write (iout,*) "Arrays UG and UGDER" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & ((ug(l,k,i),l=1,2),k=1,2), > & ((ugder(l,k,i),l=1,2),k=1,2) > enddo > write (iout,*) "Arrays UG2 and UG2DER" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & ((ug2(l,k,i),l=1,2),k=1,2), > & ((ug2der(l,k,i),l=1,2),k=1,2) > enddo > write (iout,*) "Arrays OBROT OBROT2 OBROTDER and OBROT2DER" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & (obrot(k,i),k=1,2),(obrot2(k,i),k=1,2), > & (obrot_der(k,i),k=1,2),(obrot2_der(k,i),k=1,2) > enddo > write (iout,*) "Arrays COSTAB SINTAB COSTAB2 and SINTAB2" > do i=1,nres-1 > write (iout,'(i5,4f10.5,5x,4f10.5)') i, > & costab(i),sintab(i),costab2(i),sintab2(i) > enddo > write (iout,*) "Array MUDER" > do i=1,nres-1 > write (iout,'(i5,2f10.5)') i,muder(1,i),muder(2,i) > enddo > c endif > #endif > #endif 1768a2657,2659 > #ifdef MPI > include 'mpif.h' > #endif 1770d2660 < include 'DIMENSIONS.ZSCOPT' 1771a2662 > include 'COMMON.SETUP' 1782a2674 > include 'COMMON.TIME1' 1787c2679,2681 < common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1 --- > common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, > & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, > & num_conti,j1,j2 1788a2683,2685 > #ifdef MOMENT > double precision scal_el /1.0d0/ > #else 1789a2687 > #endif 1818,1823c2716,2719 < cd if (wel_loc.gt.0.0d0) then < if (icheckgrad.eq.1) then < call vec_and_deriv_test < else < call vec_and_deriv < endif --- > c call vec_and_deriv > #ifdef TIMING > time01=MPI_Wtime() > #endif 1824a2721,2723 > #ifdef TIMING > time_mat=time_mat+MPI_Wtime()-time01 > #endif 1829c2728 < cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i) --- > cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i) 1836c2735 < num_conti_hb=0 --- > t_eelecij=0.0d0 1851a2751,2795 > c > c > c 9/27/08 AL Split the interaction loop to ensure load balancing of turn terms > C > C Loop over i,i+2 and i,i+3 pairs of the peptide groups > C > do i=iturn3_start,iturn3_end > if (itype(i).eq.21 .or. itype(i+1).eq.21 > & .or. itype(i+2).eq.21 .or. itype(i+3).eq.21) cycle > dxi=dc(1,i) > dyi=dc(2,i) > dzi=dc(3,i) > dx_normi=dc_norm(1,i) > dy_normi=dc_norm(2,i) > dz_normi=dc_norm(3,i) > xmedi=c(1,i)+0.5d0*dxi > ymedi=c(2,i)+0.5d0*dyi > zmedi=c(3,i)+0.5d0*dzi > num_conti=0 > call eelecij(i,i+2,ees,evdw1,eel_loc) > if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3) > num_cont_hb(i)=num_conti > enddo > do i=iturn4_start,iturn4_end > if (itype(i).eq.21 .or. itype(i+1).eq.21 > & .or. itype(i+3).eq.21 > & .or. itype(i+4).eq.21) cycle > dxi=dc(1,i) > dyi=dc(2,i) > dzi=dc(3,i) > dx_normi=dc_norm(1,i) > dy_normi=dc_norm(2,i) > dz_normi=dc_norm(3,i) > xmedi=c(1,i)+0.5d0*dxi > ymedi=c(2,i)+0.5d0*dyi > zmedi=c(3,i)+0.5d0*dzi > num_conti=num_cont_hb(i) > call eelecij(i,i+3,ees,evdw1,eel_loc) > if (wturn4.gt.0.0d0 .and. itype(i+2).ne.21) > & call eturn4(i,eello_turn4) > num_cont_hb(i)=num_conti > enddo ! i > c > c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 > c 1854d2797 < if (itel(i).eq.0) goto 1215 1864d2806 < num_conti=0 1865a2808 > num_conti=num_cont_hb(i) 1866a2810 > c write (iout,*) i,j,itype(i),itype(j) 1868,1869c2812,2866 < if (itel(j).eq.0) goto 1216 < ind=ind+1 --- > call eelecij(i,j,ees,evdw1,eel_loc) > enddo ! j > num_cont_hb(i)=num_conti > enddo ! i > c write (iout,*) "Number of loop steps in EELEC:",ind > cd do i=1,nres > cd write (iout,'(i3,3f10.5,5x,3f10.5)') > cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i) > cd enddo > c 12/7/99 Adam eello_turn3 will be considered as a separate energy term > ccc eel_loc=eel_loc+eello_turn3 > cd print *,"Processor",fg_rank," t_eelecij",t_eelecij > return > end > C------------------------------------------------------------------------------- > subroutine eelecij(i,j,ees,evdw1,eel_loc) > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > #ifdef MPI > include "mpif.h" > #endif > include 'COMMON.CONTROL' > include 'COMMON.IOUNITS' > include 'COMMON.GEO' > include 'COMMON.VAR' > include 'COMMON.LOCAL' > include 'COMMON.CHAIN' > include 'COMMON.DERIV' > include 'COMMON.INTERACT' > include 'COMMON.CONTACTS' > include 'COMMON.TORSION' > include 'COMMON.VECTORS' > include 'COMMON.FFIELD' > include 'COMMON.TIME1' > dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3), > & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3) > double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4), > & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4) > common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, > & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, > & num_conti,j1,j2 > c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions > #ifdef MOMENT > double precision scal_el /1.0d0/ > #else > double precision scal_el /0.5d0/ > #endif > C 12/13/98 > C 13-go grudnia roku pamietnego... > double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, > & 0.0d0,1.0d0,0.0d0, > & 0.0d0,0.0d0,1.0d0/ > c time00=MPI_Wtime() > cd write (iout,*) "eelecij",i,j > c ind=ind+1 1875,1880d2871 < C Diagnostics only!!! < c aaa=0.0D0 < c bbb=0.0D0 < c ael6i=0.0D0 < c ael3i=0.0D0 < C End diagnostics 1912d2902 < c write (iout,*) "i",i,iteli," j",j,itelj," eesij",eesij 1920a2911,2916 > > if (energy_dec) then > write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij > write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij > endif > 1925c2921 < facvdw=-6*rrmij*(ev1+evdwij) --- > facvdw=-6*rrmij*(ev1+evdwij) 1931d2926 < if (calc_grad) then 1934c2929 < * --- > * 1937a2933,2938 > c do k=1,3 > c ghalf=0.5D0*ggg(k) > c gelc(k,i)=gelc(k,i)+ghalf > c gelc(k,j)=gelc(k,j)+ghalf > c enddo > c 9/28/08 AL Gradient compotents will be summed only at the end 1939,1941c2940,2941 < ghalf=0.5D0*ggg(k) < gelc(k,i)=gelc(k,i)+ghalf < gelc(k,j)=gelc(k,j)+ghalf --- > gelc_long(k,j)=gelc_long(k,j)+ggg(k) > gelc_long(k,i)=gelc_long(k,i)-ggg(k) 1946,1950c2946,2950 < do k=i+1,j-1 < do l=1,3 < gelc(l,k)=gelc(l,k)+ggg(l) < enddo < enddo --- > cgrad do k=i+1,j-1 > cgrad do l=1,3 > cgrad gelc(l,k)=gelc(l,k)+ggg(l) > cgrad enddo > cgrad enddo 1953a2954,2959 > c do k=1,3 > c ghalf=0.5D0*ggg(k) > c gvdwpp(k,i)=gvdwpp(k,i)+ghalf > c gvdwpp(k,j)=gvdwpp(k,j)+ghalf > c enddo > c 9/28/08 AL Gradient compotents will be summed only at the end 1955,1957c2961,2962 < ghalf=0.5D0*ggg(k) < gvdwpp(k,i)=gvdwpp(k,i)+ghalf < gvdwpp(k,j)=gvdwpp(k,j)+ghalf --- > gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) > gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) 1962,1966c2967,2971 < do k=i+1,j-1 < do l=1,3 < gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) < enddo < enddo --- > cgrad do k=i+1,j-1 > cgrad do l=1,3 > cgrad gvdwpp(l,k)=gvdwpp(l,k)+ggg(l) > cgrad enddo > cgrad enddo 1975d2979 < if (calc_grad) then 1981a2986,2991 > c do k=1,3 > c ghalf=0.5D0*ggg(k) > c gelc(k,i)=gelc(k,i)+ghalf > c gelc(k,j)=gelc(k,j)+ghalf > c enddo > c 9/28/08 AL Gradient compotents will be summed only at the end 1983,1985c2993,2994 < ghalf=0.5D0*ggg(k) < gelc(k,i)=gelc(k,i)+ghalf < gelc(k,j)=gelc(k,j)+ghalf --- > gelc_long(k,j)=gelc(k,j)+ggg(k) > gelc_long(k,i)=gelc(k,i)-ggg(k) 1990,1993c2999,3010 < do k=i+1,j-1 < do l=1,3 < gelc(l,k)=gelc(l,k)+ggg(l) < enddo --- > cgrad do k=i+1,j-1 > cgrad do l=1,3 > cgrad gelc(l,k)=gelc(l,k)+ggg(l) > cgrad enddo > cgrad enddo > c 9/28/08 AL Gradient compotents will be summed only at the end > ggg(1)=facvdw*xj > ggg(2)=facvdw*yj > ggg(3)=facvdw*zj > do k=1,3 > gvdwpp(k,j)=gvdwpp(k,j)+ggg(k) > gvdwpp(k,i)=gvdwpp(k,i)-ggg(k) 2012a3030,3043 > c do k=1,3 > c ghalf=0.5D0*ggg(k) > c gelc(k,i)=gelc(k,i)+ghalf > c & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i)) > c & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1) > c gelc(k,j)=gelc(k,j)+ghalf > c & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j)) > c & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1) > c enddo > cgrad do k=i+1,j-1 > cgrad do l=1,3 > cgrad gelc(l,k)=gelc(l,k)+ggg(l) > cgrad enddo > cgrad enddo 2014,2015c3045 < ghalf=0.5D0*ggg(k) < gelc(k,i)=gelc(k,i)+ghalf --- > gelc(k,i)=gelc(k,i) 2018c3048 < gelc(k,j)=gelc(k,j)+ghalf --- > gelc(k,j)=gelc(k,j) 2020a3051,3052 > gelc_long(k,j)=gelc_long(k,j)+ggg(k) > gelc_long(k,i)=gelc_long(k,i)-ggg(k) 2022,2028d3053 < do k=i+1,j-1 < do l=1,3 < gelc(l,k)=gelc(l,k)+ggg(l) < enddo < enddo < endif < 2064,2068d3088 < C For diagnostics only < cd a22=1.0d0 < cd a23=1.0d0 < cd a32=1.0d0 < cd a33=1.0d0 2070,2072d3089 < cd write (2,*) 'fac=',fac < C For diagnostics only < cd fac=1.0d0 2080,2081c3097,3098 < cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') (uy(k,i),k=1,3), < cd & (uz(k,i),k=1,3),(uy(k,j),k=1,3),(uz(k,j),k=1,3) --- > cd write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i), > cd & uy(:,j),uz(:,j) 2086c3103 < cd write (iout,'(2i3,9f10.5/)') i,j, --- > cd write (iout,'(9f10.5/)') 2088d3104 < if (calc_grad) then 2090,2095c3106 < call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1)) < cd do k=1,3 < cd do l=1,3 < cd erder(k,l)=0.0d0 < cd enddo < cd enddo --- > call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1)) 2110,2117d3120 < cd do k=1,3 < cd do l=1,3 < cd uryg(k,l)=0.0d0 < cd urzg(k,l)=0.0d0 < cd vryg(k,l)=0.0d0 < cd vrzg(k,l)=0.0d0 < cd enddo < cd enddo 2124,2127d3126 < cd a22der=0.0d0 < cd a23der=0.0d0 < cd a32der=0.0d0 < cd a33der=0.0d0 2150,2153c3149,3152 < ghalf1=0.5d0*agg(k,1) < ghalf2=0.5d0*agg(k,2) < ghalf3=0.5d0*agg(k,3) < ghalf4=0.5d0*agg(k,4) --- > cgrad ghalf1=0.5d0*agg(k,1) > cgrad ghalf2=0.5d0*agg(k,2) > cgrad ghalf3=0.5d0*agg(k,3) > cgrad ghalf4=0.5d0*agg(k,4) 2155c3154 < & -3.0d0*uryg(k,2)*vry)+ghalf1 --- > & -3.0d0*uryg(k,2)*vry)!+ghalf1 2157c3156 < & -3.0d0*uryg(k,2)*vrz)+ghalf2 --- > & -3.0d0*uryg(k,2)*vrz)!+ghalf2 2159c3158 < & -3.0d0*urzg(k,2)*vry)+ghalf3 --- > & -3.0d0*urzg(k,2)*vry)!+ghalf3 2161c3160 < & -3.0d0*urzg(k,2)*vrz)+ghalf4 --- > & -3.0d0*urzg(k,2)*vrz)!+ghalf4 2164c3163 < & -3.0d0*uryg(k,3)*vry)+agg(k,1) --- > & -3.0d0*uryg(k,3)*vry)!+agg(k,1) 2166c3165 < & -3.0d0*uryg(k,3)*vrz)+agg(k,2) --- > & -3.0d0*uryg(k,3)*vrz)!+agg(k,2) 2168c3167 < & -3.0d0*urzg(k,3)*vry)+agg(k,3) --- > & -3.0d0*urzg(k,3)*vry)!+agg(k,3) 2170c3169 < & -3.0d0*urzg(k,3)*vrz)+agg(k,4) --- > & -3.0d0*urzg(k,3)*vrz)!+agg(k,4) 2173c3172 < & -3.0d0*vryg(k,2)*ury)+ghalf1 --- > & -3.0d0*vryg(k,2)*ury)!+ghalf1 2175c3174 < & -3.0d0*vrzg(k,2)*ury)+ghalf2 --- > & -3.0d0*vrzg(k,2)*ury)!+ghalf2 2177c3176 < & -3.0d0*vryg(k,2)*urz)+ghalf3 --- > & -3.0d0*vryg(k,2)*urz)!+ghalf3 2179c3178 < & -3.0d0*vrzg(k,2)*urz)+ghalf4 --- > & -3.0d0*vrzg(k,2)*urz)!+ghalf4 2189,2213c3188,3192 < cd aggi(k,1)=ghalf1 < cd aggi(k,2)=ghalf2 < cd aggi(k,3)=ghalf3 < cd aggi(k,4)=ghalf4 < C Derivatives in DC(i+1) < cd aggi1(k,1)=agg(k,1) < cd aggi1(k,2)=agg(k,2) < cd aggi1(k,3)=agg(k,3) < cd aggi1(k,4)=agg(k,4) < C Derivatives in DC(j) < cd aggj(k,1)=ghalf1 < cd aggj(k,2)=ghalf2 < cd aggj(k,3)=ghalf3 < cd aggj(k,4)=ghalf4 < C Derivatives in DC(j+1) < cd aggj1(k,1)=0.0d0 < cd aggj1(k,2)=0.0d0 < cd aggj1(k,3)=0.0d0 < cd aggj1(k,4)=0.0d0 < if (j.eq.nres-1 .and. i.lt.j-2) then < do l=1,4 < aggj1(k,l)=aggj1(k,l)+agg(k,l) < cd aggj1(k,l)=agg(k,l) < enddo < endif --- > cgrad if (j.eq.nres-1 .and. i.lt.j-2) then > cgrad do l=1,4 > cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l) > cgrad enddo > cgrad endif 2215,2217d3193 < endif < c goto 11111 < C Check the loc-el terms by numerical integration 2261d3236 < 11111 continue 2267c3242,3245 < cd write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3) --- > > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'eelloc',i,j,eel_loc_ij > 2270d3247 < if (calc_grad) then 2278,2284d3254 < cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij) < cd write(iout,*) 'agg ',agg < cd write(iout,*) 'aggi ',aggi < cd write(iout,*) 'aggi1',aggi1 < cd write(iout,*) 'aggj ',aggj < cd write(iout,*) 'aggj1',aggj1 < 2289,2294c3259,3269 < enddo < do k=i+2,j2 < do l=1,3 < gel_loc(l,k)=gel_loc(l,k)+ggg(l) < enddo < enddo --- > gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l) > gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l) > cgrad ghalf=0.5d0*ggg(l) > cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf > cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf > enddo > cgrad do k=i+1,j2 > cgrad do l=1,3 > cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l) > cgrad enddo > cgrad enddo 2306d3280 < endif 2308,2315d3281 < if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then < C Contributions from turns < a_temp(1,1)=a22 < a_temp(1,2)=a23 < a_temp(2,1)=a32 < a_temp(2,2)=a33 < call eturn34(i,j,eello_turn3,eello_turn4) < endif 2317c3283,3286 < if (j.gt.i+1 .and. num_conti.le.maxconts) then --- > c if (j.gt.i+1 .and. num_conti.le.maxconts) then > if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0 > & .and. num_conti.le.maxconts) then > c write (iout,*) i,j," entered corr" 2334a3304,3305 > cd write (iout,*) "i",i," j",j," num_conti",num_conti, > cd & " jcont_hb",jcont_hb(num_conti,i) 2350,2382d3320 < c if (i.eq.1) then < c a_chuj(1,1,num_conti,i)=-0.61d0 < c a_chuj(1,2,num_conti,i)= 0.4d0 < c a_chuj(2,1,num_conti,i)= 0.65d0 < c a_chuj(2,2,num_conti,i)= 0.50d0 < c else if (i.eq.2) then < c a_chuj(1,1,num_conti,i)= 0.0d0 < c a_chuj(1,2,num_conti,i)= 0.0d0 < c a_chuj(2,1,num_conti,i)= 0.0d0 < c a_chuj(2,2,num_conti,i)= 0.0d0 < c endif < C --- and its gradients < cd write (iout,*) 'i',i,' j',j < cd do kkk=1,3 < cd write (iout,*) 'iii 1 kkk',kkk < cd write (iout,*) agg(kkk,:) < cd enddo < cd do kkk=1,3 < cd write (iout,*) 'iii 2 kkk',kkk < cd write (iout,*) aggi(kkk,:) < cd enddo < cd do kkk=1,3 < cd write (iout,*) 'iii 3 kkk',kkk < cd write (iout,*) aggi1(kkk,:) < cd enddo < cd do kkk=1,3 < cd write (iout,*) 'iii 4 kkk',kkk < cd write (iout,*) aggj(kkk,:) < cd enddo < cd do kkk=1,3 < cd write (iout,*) 'iii 5 kkk',kkk < cd write (iout,*) aggj1(kkk,:) < cd enddo 2393,2395d3330 < c do mm=1,5 < c a_chuj_der(k,l,m,mm,num_conti,i)=0.0d0 < c enddo 2408,2409c3343,3356 < ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) < ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) --- > c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1) > ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1 > if (ees0tmp.gt.0) then > ees0pij=dsqrt(ees0tmp) > else > ees0pij=0 > endif > c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2) > ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2 > if (ees0tmp.gt.0) then > ees0mij=dsqrt(ees0tmp) > else > ees0mij=0 > endif 2418,2421c3365,3366 < c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij, < c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont < facont_hb(num_conti,i)=fcont < if (calc_grad) then --- > c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij, > c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont 2448a3394 > facont_hb(num_conti,i)=fcont 2472,2474c3418,3424 < ghalfp=0.5D0*gggp(k) < ghalfm=0.5D0*gggm(k) < gacontp_hb1(k,num_conti,i)=ghalfp --- > c > c 10/24/08 cgrad and ! comments indicate the parts of the code removed > c following the change of gradient-summation algorithm. > c > cgrad ghalfp=0.5D0*gggp(k) > cgrad ghalfm=0.5D0*gggm(k) > gacontp_hb1(k,num_conti,i)=!ghalfp 2477c3427 < gacontp_hb2(k,num_conti,i)=ghalfp --- > gacontp_hb2(k,num_conti,i)=!ghalfp 2481c3431 < gacontm_hb1(k,num_conti,i)=ghalfm --- > gacontm_hb1(k,num_conti,i)=!ghalfm 2484c3434 < gacontm_hb2(k,num_conti,i)=ghalfm --- > gacontm_hb2(k,num_conti,i)=!ghalfm 2489d3438 < endif 2503,2513c3452,3469 < 1216 continue < enddo ! j < num_cont_hb(i)=num_conti < 1215 continue < enddo ! i < cd do i=1,nres < cd write (iout,'(i3,3f10.5,5x,3f10.5)') < cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i) < cd enddo < c 12/7/99 Adam eello_turn3 will be considered as a separate energy term < ccc eel_loc=eel_loc+eello_turn3 --- > if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then > do k=1,4 > do l=1,3 > ghalf=0.5d0*agg(l,k) > aggi(l,k)=aggi(l,k)+ghalf > aggi1(l,k)=aggi1(l,k)+agg(l,k) > aggj(l,k)=aggj(l,k)+ghalf > enddo > enddo > if (j.eq.nres-1 .and. i.lt.j-2) then > do k=1,4 > do l=1,3 > aggj1(l,k)=aggj1(l,k)+agg(l,k) > enddo > enddo > endif > endif > c t_eelecij=t_eelecij+MPI_Wtime()-time00 2517c3473 < subroutine eturn34(i,j,eello_turn3,eello_turn4) --- > subroutine eturn3(i,eello_turn3) 2521d3476 < include 'DIMENSIONS.ZSCOPT' 2532a3488 > include 'COMMON.CONTROL' 2538,2540c3494,3503 < & aggj(3,4),aggj1(3,4),a_temp(2,2) < common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2 < if (j.eq.i+2) then --- > & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) > common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, > & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, > & num_conti,j1,j2 > j=i+2 > c write (iout,*) "eturn3",i,j,j1,j2 > a_temp(1,1)=a22 > a_temp(1,2)=a23 > a_temp(2,1)=a32 > a_temp(2,2)=a33 2555a3519,3520 > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) 2559d3523 < if (calc_grad) then 2562,2563c3526,3527 < call transpose2(auxmat2(1,1),pizda(1,1)) < call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) --- > call transpose2(auxmat2(1,1),auxmat3(1,1)) > call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) 2567,2568c3531,3532 < call transpose2(auxmat2(1,1),pizda(1,1)) < call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) --- > call transpose2(auxmat2(1,1),auxmat3(1,1)) > call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1)) 2573,2576c3537,3544 < a_temp(1,1)=aggi(l,1) < a_temp(1,2)=aggi(l,2) < a_temp(2,1)=aggi(l,3) < a_temp(2,2)=aggi(l,4) --- > c ghalf1=0.5d0*agg(l,1) > c ghalf2=0.5d0*agg(l,2) > c ghalf3=0.5d0*agg(l,3) > c ghalf4=0.5d0*agg(l,4) > a_temp(1,1)=aggi(l,1)!+ghalf1 > a_temp(1,2)=aggi(l,2)!+ghalf2 > a_temp(2,1)=aggi(l,3)!+ghalf3 > a_temp(2,2)=aggi(l,4)!+ghalf4 2580,2583c3548,3551 < a_temp(1,1)=aggi1(l,1) < a_temp(1,2)=aggi1(l,2) < a_temp(2,1)=aggi1(l,3) < a_temp(2,2)=aggi1(l,4) --- > a_temp(1,1)=aggi1(l,1)!+agg(l,1) > a_temp(1,2)=aggi1(l,2)!+agg(l,2) > a_temp(2,1)=aggi1(l,3)!+agg(l,3) > a_temp(2,2)=aggi1(l,4)!+agg(l,4) 2587,2590c3555,3558 < a_temp(1,1)=aggj(l,1) < a_temp(1,2)=aggj(l,2) < a_temp(2,1)=aggj(l,3) < a_temp(2,2)=aggj(l,4) --- > a_temp(1,1)=aggj(l,1)!+ghalf1 > a_temp(1,2)=aggj(l,2)!+ghalf2 > a_temp(2,1)=aggj(l,3)!+ghalf3 > a_temp(2,2)=aggj(l,4)!+ghalf4 2602,2603c3570,3598 < endif < else if (j.eq.i+3 .and. itype(i+2).ne.21) then --- > return > end > C------------------------------------------------------------------------------- > subroutine eturn4(i,eello_turn4) > C Third- and fourth-order contributions from turns > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'COMMON.IOUNITS' > include 'COMMON.GEO' > include 'COMMON.VAR' > include 'COMMON.LOCAL' > include 'COMMON.CHAIN' > include 'COMMON.DERIV' > include 'COMMON.INTERACT' > include 'COMMON.CONTACTS' > include 'COMMON.TORSION' > include 'COMMON.VECTORS' > include 'COMMON.FFIELD' > include 'COMMON.CONTROL' > dimension ggg(3) > double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2), > & e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2), > & e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2) > double precision agg(3,4),aggi(3,4),aggi1(3,4), > & aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2) > common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33, > & dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi, > & num_conti,j1,j2 > j=i+3 2615a3611,3615 > c write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2 > a_temp(1,1)=a22 > a_temp(1,2)=a23 > a_temp(2,1)=a32 > a_temp(2,2)=a33 2618a3619 > c write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 2631a3633,3634 > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'eturn4',i,j,-(s1+s2+s3) 2635d3637 < if (calc_grad) then 2658,2659c3660,3661 < call matmat2(auxmat(1,1),e2t(1,1),auxmat(1,1)) < call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1)) --- > call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1)) > call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1)) 2739a3742 > c write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3 2742,2743d3744 < endif < endif 2779a3781,3877 > subroutine escp_soft_sphere(evdw2,evdw2_14) > C > C This subroutine calculates the excluded-volume interaction energy between > C peptide-group centers and side chains and its gradient in virtual-bond and > C side-chain vectors. > C > implicit real*8 (a-h,o-z) > include 'DIMENSIONS' > include 'COMMON.GEO' > include 'COMMON.VAR' > include 'COMMON.LOCAL' > include 'COMMON.CHAIN' > include 'COMMON.DERIV' > include 'COMMON.INTERACT' > include 'COMMON.FFIELD' > include 'COMMON.IOUNITS' > include 'COMMON.CONTROL' > dimension ggg(3) > evdw2=0.0D0 > evdw2_14=0.0d0 > r0_scp=4.5d0 > cd print '(a)','Enter ESCP' > cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e > do i=iatscp_s,iatscp_e > if (itype(i).eq.21 .or. itype(i+1).eq.21) cycle > iteli=itel(i) > xi=0.5D0*(c(1,i)+c(1,i+1)) > yi=0.5D0*(c(2,i)+c(2,i+1)) > zi=0.5D0*(c(3,i)+c(3,i+1)) > > do iint=1,nscp_gr(i) > > do j=iscpstart(i,iint),iscpend(i,iint) > if (itype(j).eq.21) cycle > itypj=itype(j) > C Uncomment following three lines for SC-p interactions > c xj=c(1,nres+j)-xi > c yj=c(2,nres+j)-yi > c zj=c(3,nres+j)-zi > C Uncomment following three lines for Ca-p interactions > xj=c(1,j)-xi > yj=c(2,j)-yi > zj=c(3,j)-zi > rij=xj*xj+yj*yj+zj*zj > r0ij=r0_scp > r0ijsq=r0ij*r0ij > if (rij.lt.r0ijsq) then > evdwij=0.25d0*(rij-r0ijsq)**2 > fac=rij-r0ijsq > else > evdwij=0.0d0 > fac=0.0d0 > endif > evdw2=evdw2+evdwij > C > C Calculate contributions to the gradient in the virtual-bond and SC vectors. > C > ggg(1)=xj*fac > ggg(2)=yj*fac > ggg(3)=zj*fac > cgrad if (j.lt.i) then > cd write (iout,*) 'j C Uncomment following three lines for SC-p interactions > c do k=1,3 > c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) > c enddo > cgrad else > cd write (iout,*) 'j>i' > cgrad do k=1,3 > cgrad ggg(k)=-ggg(k) > C Uncomment following line for SC-p interactions > c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k) > cgrad enddo > cgrad endif > cgrad do k=1,3 > cgrad gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k) > cgrad enddo > cgrad kstart=min0(i+1,j) > cgrad kend=max0(i-1,j-1) > cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend > cd write (iout,*) ggg(1),ggg(2),ggg(3) > cgrad do k=kstart,kend > cgrad do l=1,3 > cgrad gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l) > cgrad enddo > cgrad enddo > do k=1,3 > gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) > gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) > enddo > enddo > > enddo ! iint > enddo ! i > return > end > C----------------------------------------------------------------------------- 2788d3885 < include 'DIMENSIONS.ZSCOPT' 2796a3894 > include 'COMMON.CONTROL' 2801,2802c3899 < c write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e, < c & ' scal14',scal14 --- > cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e 2806,2808d3902 < c write (iout,*) "i",i," iteli",iteli," nscp_gr",nscp_gr(i), < c & " iscp",(iscpstart(i,j),iscpend(i,j),j=1,nscp_gr(i)) < if (iteli.eq.0) goto 1225 2836d3929 < c write (iout,*) i,j,evdwij 2838c3931,3932 < if (calc_grad) then --- > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'evdw2',i,j,evdwij 2846c3940 < if (j.lt.i) then --- > cgrad if (j.lt.i) then 2852c3946 < else --- > cgrad else 2854,2855c3948,3949 < do k=1,3 < ggg(k)=-ggg(k) --- > cgrad do k=1,3 > cgrad ggg(k)=-ggg(k) 2857,2864c3951,3959 < c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k) < enddo < endif < do k=1,3 < gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k) < enddo < kstart=min0(i+1,j) < kend=max0(i-1,j-1) --- > ccgrad gradx_scp(k,j)=gradx_scp(k,j)-ggg(k) > c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k) > cgrad enddo > cgrad endif > cgrad do k=1,3 > cgrad gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k) > cgrad enddo > cgrad kstart=min0(i+1,j) > cgrad kend=max0(i-1,j-1) 2867,2870c3962,3969 < do k=kstart,kend < do l=1,3 < gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l) < enddo --- > cgrad do k=kstart,kend > cgrad do l=1,3 > cgrad gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l) > cgrad enddo > cgrad enddo > do k=1,3 > gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k) > gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k) 2872d3970 < endif 2873a3972 > 2875d3973 < 1225 continue 2879a3978 > gvdwc_scpp(j,i)=expon*gvdwc_scpp(j,i) 2901d3999 < include 'DIMENSIONS.ZSCOPT' 2906a4005 > include 'COMMON.IOUNITS' 2909,2910c4008,4009 < cd print *,'edis: nhpb=',nhpb,' fbr=',fbr < cd print *,'link_start=',link_start,' link_end=',link_end --- > cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr > cd write(iout,*)'link_start=',link_start,' link_end=',link_end 2924a4024 > cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj 2929a4030 > cd write (iout,*) "eij",eij 2957,2960c4058,4065 < do j=iii,jjj-1 < do k=1,3 < ghpbc(k,j)=ghpbc(k,j)+ggg(k) < enddo --- > cgrad do j=iii,jjj-1 > cgrad do k=1,3 > cgrad ghpbc(k,j)=ghpbc(k,j)+ggg(k) > cgrad enddo > cgrad enddo > do k=1,3 > ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) > ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) 2978d4082 < include 'DIMENSIONS.ZSCOPT' 2994c4098,4099 < dsci_inv=dsc_inv(itypi) --- > c dsci_inv=dsc_inv(itypi) > dsci_inv=vbld_inv(nres+i) 2996c4101,4102 < dscj_inv=dsc_inv(itypj) --- > c dscj_inv=dsc_inv(itypj) > dscj_inv=vbld_inv(nres+j) 3034,3040c4140,4148 < gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k) < enddo < do k=1,3 < ghpbx(k,i)=ghpbx(k,i)-gg(k) < & +(eom12*dc_norm(k,nres+j)+eom1*erij(k))*dsci_inv < ghpbx(k,j)=ghpbx(k,j)+gg(k) < & +(eom12*dc_norm(k,nres+i)+eom2*erij(k))*dscj_inv --- > ggk=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k) > ghpbx(k,i)=ghpbx(k,i)-ggk > & +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) > & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv > ghpbx(k,j)=ghpbx(k,j)+ggk > & +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) > & +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv > ghpbc(k,i)=ghpbc(k,i)-ggk > ghpbc(k,j)=ghpbc(k,j)+ggk 3045,3049c4153,4157 < do k=i,j-1 < do l=1,3 < ghpbc(l,k)=ghpbc(l,k)+gg(l) < enddo < enddo --- > cgrad do k=i,j-1 > cgrad do l=1,3 > cgrad ghpbc(l,k)=ghpbc(l,k)+gg(l) > cgrad enddo > cgrad enddo 3059d4166 < include 'DIMENSIONS.ZSCOPT' 3070c4177 < logical energy_dec /.false./ --- > include 'COMMON.SETUP' 3073,3074c4180,4181 < write (iout,*) "distchainmax",distchainmax < do i=nnt+1,nct --- > estr1=0.0d0 > do i=ibondp_start,ibondp_end 3081c4188 < if (energy_dec) write(iout,*) --- > if (energy_dec) write(iout,*) 3085,3090c4192,4199 < diff = vbld(i)-vbldp0 < c write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff < estr=estr+diff*diff < do j=1,3 < gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) < enddo --- > diff = vbld(i)-vbldp0 > if (energy_dec) write (iout,*) > & "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff > estr=estr+diff*diff > do j=1,3 > gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i) > enddo > c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) 3092d4200 < 3094c4202 < estr=0.5d0*AKP*estr --- > estr=0.5d0*AKP*estr+estr1 3098c4206 < do i=nnt,nct --- > do i=ibond_start,ibond_end 3104,3105c4212,4214 < c write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff, < c & AKSC(1,iti),AKSC(1,iti)*diff*diff --- > if (energy_dec) write (iout,*) > & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff, > & AKSC(1,iti),AKSC(1,iti)*diff*diff 3112c4221 < diff=vbld(i+nres)-vbldsc0(j,iti) --- > diff=vbld(i+nres)-vbldsc0(j,iti) 3132c4241 < usumsqder=usumsqder+ud(j)*uprod2 --- > usumsqder=usumsqder+ud(j)*uprod2 3134,3135d4242 < c write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti), < c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) 3144c4251 < end --- > end 3154d4260 < include 'DIMENSIONS.ZSCOPT' 3163a4270 > include 'COMMON.CONTROL' 3169,3170c4276,4277 < time11=dexp(-2*time) < time12=1.0d0 --- > c time11=dexp(-2*time) > c time12=1.0d0 3172d4278 < c write (iout,*) "nres",nres 3174d4279 < c write (iout,*) ithet_start,ithet_end 3182,3185c4287,4288 < phii=phi(i) < icrc=0 < call proc_proc(phii,icrc) < if (icrc.eq.1) phii=150.0 --- > phii=phi(i) > if (phii.ne.phii) phii=150.0 3191c4294 < else --- > else 3197,3200c4300,4301 < phii1=phi(i+1) < icrc=0 < call proc_proc(phii1,icrc) < if (icrc.eq.1) phii1=150.0 --- > phii1=phi(i+1) > if (phii1.ne.phii1) phii1=150.0 3211c4312 < endif --- > endif 3221d4321 < c write (iout,*) "thet_pred_mean",thet_pred_mean 3224d4323 < c write (iout,*) "thet_pred_mean",thet_pred_mean 3250,3251c4349,4350 < c write (iout,'(2i3,3f8.3,f10.5)') i,it,rad2deg*theta(i), < c & rad2deg*phii,rad2deg*phii1,ethetai --- > if (energy_dec) write (iout,'(a6,i5,0pf7.3)') > & 'ebend',i,ethetai 3255d4353 < 1215 continue 3379d4476 < include 'DIMENSIONS.ZSCOPT' 3396d4492 < c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) 3449,3451d4544 < c write (iout,*) "i",i," ityp1",itype(i-2),ityp1, < c & " ityp2",itype(i-1),ityp2," ityp3",itype(i),ityp3 < c call flush(iout) 3571d4663 < include 'DIMENSIONS.ZSCOPT' 3580a4673 > include 'COMMON.CONTROL' 3598d4690 < c write (iout,*) "i",i," x",x(1),x(2),x(3) 3672c4764,4766 < c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc --- > if (energy_dec) write (iout,'(a6,i5,0pf7.3)') > & 'escloc',i,escloci > c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc 3745a4840,4844 > #ifdef OSF > adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin > if(adexp.ne.adexp) adexp=1.0 > expfac=dexp(adexp) > #else 3746a4846 > #endif 3856d4955 < include 'DIMENSIONS.ZSCOPT' 3991,3992c5090 < c write (2,*) "escloc",escloc < if (.not. calc_grad) goto 1 --- > c write (2,*) "i",i," escloc",sumene,escloc 4168a5267,5304 > c------------------------------------------------------------------------------ > double precision function enesc(x,xx,yy,zz,cost2,sint2) > implicit none > double precision x(65),xx,yy,zz,cost2,sint2,sumene1,sumene2, > & sumene3,sumene4,sumene,dsc_i,dp2_i,dscp1,dscp2,s1,s1_6,s2,s2_6 > sumene1= x(1)+ x(2)*xx+ x(3)*yy+ x(4)*zz+ x(5)*xx**2 > & + x(6)*yy**2+ x(7)*zz**2+ x(8)*xx*zz+ x(9)*xx*yy > & + x(10)*yy*zz > sumene2= x(11) + x(12)*xx + x(13)*yy + x(14)*zz + x(15)*xx**2 > & + x(16)*yy**2 + x(17)*zz**2 + x(18)*xx*zz + x(19)*xx*yy > & + x(20)*yy*zz > sumene3= x(21) +x(22)*xx +x(23)*yy +x(24)*zz +x(25)*xx**2 > & +x(26)*yy**2 +x(27)*zz**2 +x(28)*xx*zz +x(29)*xx*yy > & +x(30)*yy*zz +x(31)*xx**3 +x(32)*yy**3 +x(33)*zz**3 > & +x(34)*(xx**2)*yy +x(35)*(xx**2)*zz +x(36)*(yy**2)*xx > & +x(37)*(yy**2)*zz +x(38)*(zz**2)*xx +x(39)*(zz**2)*yy > & +x(40)*xx*yy*zz > sumene4= x(41) +x(42)*xx +x(43)*yy +x(44)*zz +x(45)*xx**2 > & +x(46)*yy**2 +x(47)*zz**2 +x(48)*xx*zz +x(49)*xx*yy > & +x(50)*yy*zz +x(51)*xx**3 +x(52)*yy**3 +x(53)*zz**3 > & +x(54)*(xx**2)*yy +x(55)*(xx**2)*zz +x(56)*(yy**2)*xx > & +x(57)*(yy**2)*zz +x(58)*(zz**2)*xx +x(59)*(zz**2)*yy > & +x(60)*xx*yy*zz > dsc_i = 0.743d0+x(61) > dp2_i = 1.9d0+x(62) > dscp1=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i > & *(xx*cost2+yy*sint2)) > dscp2=dsqrt(dsc_i**2+dp2_i**2-2*dsc_i*dp2_i > & *(xx*cost2-yy*sint2)) > s1=(1+x(63))/(0.1d0 + dscp1) > s1_6=(1+x(64))/(0.1d0 + dscp1**6) > s2=(1+x(65))/(0.1d0 + dscp2) > s2_6=(1+x(65))/(0.1d0 + dscp2**6) > sumene = ( sumene3*sint2 + sumene1)*(s1+s1_6) > & + (sumene4*cost2 +sumene2)*(s2+s2_6) > enesc=sumene > return > end 4207d5342 < include 'DIMENSIONS.ZSCOPT' 4252c5387 < subroutine etor(etors,edihcnstr,fact) --- > subroutine etor(etors,edihcnstr) 4255d5389 < include 'DIMENSIONS.ZSCOPT' 4266a5401 > include 'COMMON.CONTROL' 4273c5408,5409 < if (itype(i-2).eq.21 .or. itype(i-1).eq.21 --- > etors_ii=0.0D0 > if (itype(i-2).eq.21 .or. itype(i-1).eq.21 4286a5423 > if (energy_dec) etors_ii=etors_ii+etorsi-v1(1,3,3) 4294a5432,5433 > if (energy_dec) etors_ii=etors_ii+ > & v2ij*sinphi+dabs(v1ij)+dabs(v2ij) 4303a5443,5444 > if (energy_dec) etors_ii=etors_ii+ > & v1ij*cosphi+v2ij*sinphi+dabs(v1ij)+dabs(v2ij) 4306a5448,5449 > if (energy_dec) write (iout,'(a6,i5,0pf7.3)') > 'etor',i,etors_ii 4311c5454 < gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci --- > gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci 4335a5479,5483 > subroutine etor_d(etors_d) > etors_d=0.0d0 > return > end > c---------------------------------------------------------------------------- 4337c5485 < subroutine etor(etors,edihcnstr,fact) --- > subroutine etor(etors,edihcnstr) 4340d5487 < include 'DIMENSIONS.ZSCOPT' 4351a5499 > include 'COMMON.CONTROL' 4355c5503 < c lprn=.true. --- > c lprn=.true. 4358c5506 < if (itype(i-2).eq.21 .or. itype(i-1).eq.21 --- > if (itype(i-2).eq.21 .or. itype(i-1).eq.21 4360c5508 < if (itel(i-2).eq.0 .or. itel(i-1).eq.0) goto 1215 --- > etors_ii=0.0D0 4371a5520,5521 > if (energy_dec) etors_ii=etors_ii+ > & v1ij*cosphi+v2ij*sinphi 4387a5538,5539 > if (energy_dec) etors_ii=etors_ii+ > & vl1ij*pom1 4392a5545,5546 > if (energy_dec) write (iout,'(a6,i5,0pf7.3)') > & 'etor',i,etors_ii-v0(itori,itori1) 4397c5551 < gloc(i-3,icg)=gloc(i-3,icg)+wtor*fact*gloci --- > gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci 4399d5552 < 1215 continue 4403c5556,5557 < do i=1,ndih_constr --- > c do i=1,ndih_constr > do i=idihconstr_start,idihconstr_end 4407d5560 < edihi=0.0d0 4412d5564 < edihi=0.25d0*ftors*difi**4 4417d5568 < edihi=0.25d0*ftors*difi**4 4419c5570 < difi=0.0d0 --- > difi=0.0 4421,4424c5572,5574 < c write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi, < c & drange(i),edihi < ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, < ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) --- > cd write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, > cd & rad2deg*phi0(i), rad2deg*drange(i), > cd & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) 4426c5576 < ! write (iout,*) 'edihcnstr',edihcnstr --- > cd write (iout,*) 'edihcnstr',edihcnstr 4430c5580 < subroutine etor_d(etors_d,fact2) --- > subroutine etor_d(etors_d) 4434d5583 < include 'DIMENSIONS.ZSCOPT' 4451c5600 < do i=iphi_start,iphi_end-1 --- > do i=iphid_start,iphid_end 4454,4455d5602 < if (itel(i-2).eq.0 .or. itel(i-1).eq.0 .or. itel(i).eq.0) < & goto 1215 4496,4498c5643,5644 < gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*fact2*gloci1 < gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*fact2*gloci2 < 1215 continue --- > gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1 > gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2 4509c5655 < c of residues computed from AM1 energy surfaces of terminally-blocked --- > c of residues computed from AM1 energy surfaces of terminally-blocked 4513d5658 < include 'DIMENSIONS.ZSCOPT' 4551c5696 < gsccor_loc(i-3)=gloci --- > gsccor_loc(i-3)=gsccor_loc(i-3)+gloci 4555c5700 < c------------------------------------------------------------------------------ --- > c---------------------------------------------------------------------------- 4638,4702c5783,5794 < enddo < do m=i,j-1 < do ll=1,3 < gradcorr(ll,m)=gradcorr(ll,m)+gx(ll) < enddo < enddo < do m=k,l-1 < do ll=1,3 < gradcorr(ll,m)=gradcorr(ll,m)+gx1(ll) < enddo < enddo < esccorr=-eij*ekl < return < end < c------------------------------------------------------------------------------ < #ifdef MPL < subroutine pack_buffer(dimen1,dimen2,atom,indx,buffer) < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < integer dimen1,dimen2,atom,indx < double precision buffer(dimen1,dimen2) < double precision zapas < common /contacts_hb/ zapas(3,20,maxres,7), < & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), < & num_cont_hb(maxres),jcont_hb(20,maxres) < num_kont=num_cont_hb(atom) < do i=1,num_kont < do k=1,7 < do j=1,3 < buffer(i,indx+(k-1)*3+j)=zapas(j,i,atom,k) < enddo ! j < enddo ! k < buffer(i,indx+22)=facont_hb(i,atom) < buffer(i,indx+23)=ees0p(i,atom) < buffer(i,indx+24)=ees0m(i,atom) < buffer(i,indx+25)=dfloat(jcont_hb(i,atom)) < enddo ! i < buffer(1,indx+26)=dfloat(num_kont) < return < end < c------------------------------------------------------------------------------ < subroutine unpack_buffer(dimen1,dimen2,atom,indx,buffer) < implicit real*8 (a-h,o-z) < include 'DIMENSIONS' < integer dimen1,dimen2,atom,indx < double precision buffer(dimen1,dimen2) < double precision zapas < common /contacts_hb/ zapas(3,20,maxres,7), < & facont_hb(20,maxres),ees0p(20,maxres),ees0m(20,maxres), < & num_cont_hb(maxres),jcont_hb(20,maxres) < num_kont=buffer(1,indx+26) < 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,7 < do j=1,3 < zapas(j,ii,atom,k)=buffer(i,indx+(k-1)*3+j) < enddo ! j < enddo ! k < facont_hb(ii,atom)=buffer(i,indx+22) < ees0p(ii,atom)=buffer(i,indx+23) < ees0m(ii,atom)=buffer(i,indx+24) < jcont_hb(ii,atom)=buffer(i,indx+25) < enddo ! i --- > enddo > do m=i,j-1 > do ll=1,3 > gradcorr(ll,m)=gradcorr(ll,m)+gx(ll) > enddo > enddo > do m=k,l-1 > do ll=1,3 > gradcorr(ll,m)=gradcorr(ll,m)+gx1(ll) > enddo > enddo > esccorr=-eij*ekl 4706d5797 < #endif 4711d5801 < include 'DIMENSIONS.ZSCOPT' 4713,4714c5803,5812 < #ifdef MPL < include 'COMMON.INFO' --- > #ifdef MPI > include "mpif.h" > parameter (max_cont=maxconts) > parameter (max_dim=26) > integer source,CorrelType,CorrelID,CorrelType1,CorrelID1,Error > double precision zapas(max_dim,maxconts,max_fg_procs), > & zapas_recv(max_dim,maxconts,max_fg_procs) > common /przechowalnia/ zapas > integer status(MPI_STATUS_SIZE),req(maxconts*2), > & status_array(MPI_STATUS_SIZE,maxconts*2) 4715a5814 > include 'COMMON.SETUP' 4720,4728c5819,5821 < #ifdef MPL < parameter (max_cont=maxconts) < parameter (max_dim=2*(8*3+2)) < parameter (msglen1=max_cont*max_dim*4) < parameter (msglen2=2*msglen1) < integer source,CorrelType,CorrelID,Error < double precision buffer(max_cont,max_dim) < #endif < double precision gx(3),gx1(3) --- > include 'COMMON.CONTROL' > include 'COMMON.LOCAL' > double precision gx(3),gx1(3),time00 4733c5826 < #ifdef MPL --- > #ifdef MPI 4736c5829 < if (fgProcs.le.1) goto 30 --- > if (nfgtasks.le.1) goto 30 4738c5831 < write (iout,'(a)') 'Contact function values:' --- > write (iout,'(a)') 'Contact function values before RECEIVE:' 4745,4746c5838,5917 < C Caution! Following code assumes that electrostatic interactions concerning < C a given atom are split among at most two processors! --- > call flush(iout) > do i=1,ntask_cont_from > ncont_recv(i)=0 > enddo > do i=1,ntask_cont_to > ncont_sent(i)=0 > enddo > c write (iout,*) "ntask_cont_from",ntask_cont_from," ntask_cont_to", > c & ntask_cont_to > C Make the list of contacts to send to send to other procesors > c write (iout,*) "limits",max0(iturn4_end-1,iatel_s),iturn3_end > c call flush(iout) > do i=iturn3_start,iturn3_end > c write (iout,*) "make contact list turn3",i," num_cont", > c & num_cont_hb(i) > call add_hb_contact(i,i+2,iturn3_sent_local(1,i)) > enddo > do i=iturn4_start,iturn4_end > c write (iout,*) "make contact list turn4",i," num_cont", > c & num_cont_hb(i) > call add_hb_contact(i,i+3,iturn4_sent_local(1,i)) > enddo > do ii=1,nat_sent > i=iat_sent(ii) > c write (iout,*) "make contact list longrange",i,ii," num_cont", > c & num_cont_hb(i) > do j=1,num_cont_hb(i) > do k=1,4 > jjc=jcont_hb(j,i) > iproc=iint_sent_local(k,jjc,ii) > c write (iout,*) "i",i," j",j," k",k," jjc",jjc," iproc",iproc > if (iproc.gt.0) then > ncont_sent(iproc)=ncont_sent(iproc)+1 > nn=ncont_sent(iproc) > zapas(1,nn,iproc)=i > zapas(2,nn,iproc)=jjc > zapas(3,nn,iproc)=facont_hb(j,i) > zapas(4,nn,iproc)=ees0p(j,i) > zapas(5,nn,iproc)=ees0m(j,i) > zapas(6,nn,iproc)=gacont_hbr(1,j,i) > zapas(7,nn,iproc)=gacont_hbr(2,j,i) > zapas(8,nn,iproc)=gacont_hbr(3,j,i) > zapas(9,nn,iproc)=gacontm_hb1(1,j,i) > zapas(10,nn,iproc)=gacontm_hb1(2,j,i) > zapas(11,nn,iproc)=gacontm_hb1(3,j,i) > zapas(12,nn,iproc)=gacontp_hb1(1,j,i) > zapas(13,nn,iproc)=gacontp_hb1(2,j,i) > zapas(14,nn,iproc)=gacontp_hb1(3,j,i) > zapas(15,nn,iproc)=gacontm_hb2(1,j,i) > zapas(16,nn,iproc)=gacontm_hb2(2,j,i) > zapas(17,nn,iproc)=gacontm_hb2(3,j,i) > zapas(18,nn,iproc)=gacontp_hb2(1,j,i) > zapas(19,nn,iproc)=gacontp_hb2(2,j,i) > zapas(20,nn,iproc)=gacontp_hb2(3,j,i) > zapas(21,nn,iproc)=gacontm_hb3(1,j,i) > zapas(22,nn,iproc)=gacontm_hb3(2,j,i) > zapas(23,nn,iproc)=gacontm_hb3(3,j,i) > zapas(24,nn,iproc)=gacontp_hb3(1,j,i) > zapas(25,nn,iproc)=gacontp_hb3(2,j,i) > zapas(26,nn,iproc)=gacontp_hb3(3,j,i) > endif > enddo > enddo > enddo > if (lprn) then > write (iout,*) > & "Numbers of contacts to be sent to other processors", > & (ncont_sent(i),i=1,ntask_cont_to) > write (iout,*) "Contacts sent" > do ii=1,ntask_cont_to > nn=ncont_sent(ii) > iproc=itask_cont_to(ii) > write (iout,*) nn," contacts to processor",iproc, > & " of CONT_TO_COMM group" > do i=1,nn > write(iout,'(2f5.0,4f10.5)')(zapas(j,i,ii),j=1,5) > enddo > enddo > call flush(iout) > endif 4748,4836c5919,6040 < CorrelID=MyID+1 < ldone=.false. < do i=1,max_cont < do j=1,max_dim < buffer(i,j)=0.0D0 < enddo < enddo < mm=mod(MyRank,2) < cd write (iout,*) 'MyRank',MyRank,' mm',mm < if (mm) 20,20,10 < 10 continue < cd write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone < if (MyRank.gt.0) then < C Send correlation contributions to the preceding processor < msglen=msglen1 < nn=num_cont_hb(iatel_s) < call pack_buffer(max_cont,max_dim,iatel_s,0,buffer) < cd write (iout,*) 'The BUFFER array:' < cd do i=1,nn < cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26) < cd enddo < if (ielstart(iatel_s).gt.iatel_s+ispp) then < msglen=msglen2 < call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer) < C Clear the contacts of the atom passed to the neighboring processor < nn=num_cont_hb(iatel_s+1) < cd do i=1,nn < cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26) < cd enddo < num_cont_hb(iatel_s)=0 < endif < cd write (iout,*) 'Processor ',MyID,MyRank, < cd & ' is sending correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen < cd write (*,*) 'Processor ',MyID,MyRank, < cd & ' is sending correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen,' CorrelType=',CorrelType < call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID) < cd write (iout,*) 'Processor ',MyID, < cd & ' has sent correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen,' CorrelID=',CorrelID < cd write (*,*) 'Processor ',MyID, < cd & ' has sent correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen,' CorrelID=',CorrelID < msglen=msglen1 < endif ! (MyRank.gt.0) < if (ldone) goto 30 < ldone=.true. < 20 continue < cd write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone < if (MyRank.lt.fgProcs-1) then < C Receive correlation contributions from the next processor < msglen=msglen1 < if (ielend(iatel_e).lt.nct-1) msglen=msglen2 < cd write (iout,*) 'Processor',MyID, < cd & ' is receiving correlation contribution from processor',MyID+1, < cd & ' msglen=',msglen,' CorrelType=',CorrelType < cd write (*,*) 'Processor',MyID, < cd & ' is receiving correlation contribution from processor',MyID+1, < cd & ' msglen=',msglen,' CorrelType=',CorrelType < nbytes=-1 < do while (nbytes.le.0) < call mp_probe(MyID+1,CorrelType,nbytes) < enddo < cd print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes < call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes) < cd write (iout,*) 'Processor',MyID, < cd & ' has received correlation contribution from processor',MyID+1, < cd & ' msglen=',msglen,' nbytes=',nbytes < cd write (iout,*) 'The received BUFFER array:' < cd do i=1,max_cont < cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52) < cd enddo < if (msglen.eq.msglen1) then < call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer) < else if (msglen.eq.msglen2) then < call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer) < call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer) < else < write (iout,*) < & 'ERROR!!!! message length changed while processing correlations.' < write (*,*) < & 'ERROR!!!! message length changed while processing correlations.' < call mp_stopall(Error) < endif ! msglen.eq.msglen1 < endif ! MyRank.lt.fgProcs-1 < if (ldone) goto 30 < ldone=.true. < goto 10 --- > CorrelID=fg_rank+1 > CorrelType1=478 > CorrelID1=nfgtasks+fg_rank+1 > ireq=0 > C Receive the numbers of needed contacts from other processors > do ii=1,ntask_cont_from > iproc=itask_cont_from(ii) > ireq=ireq+1 > call MPI_Irecv(ncont_recv(ii),1,MPI_INTEGER,iproc,CorrelType, > & FG_COMM,req(ireq),IERR) > enddo > c write (iout,*) "IRECV ended" > c call flush(iout) > C Send the number of contacts needed by other processors > do ii=1,ntask_cont_to > iproc=itask_cont_to(ii) > ireq=ireq+1 > call MPI_Isend(ncont_sent(ii),1,MPI_INTEGER,iproc,CorrelType, > & FG_COMM,req(ireq),IERR) > enddo > c write (iout,*) "ISEND ended" > c write (iout,*) "number of requests (nn)",ireq > call flush(iout) > if (ireq.gt.0) > & call MPI_Waitall(ireq,req,status_array,ierr) > c write (iout,*) > c & "Numbers of contacts to be received from other processors", > c & (ncont_recv(i),i=1,ntask_cont_from) > c call flush(iout) > C Receive contacts > ireq=0 > do ii=1,ntask_cont_from > iproc=itask_cont_from(ii) > nn=ncont_recv(ii) > c write (iout,*) "Receiving",nn," contacts from processor",iproc, > c & " of CONT_TO_COMM group" > call flush(iout) > if (nn.gt.0) then > ireq=ireq+1 > call MPI_Irecv(zapas_recv(1,1,ii),nn*max_dim, > & MPI_DOUBLE_PRECISION,iproc,CorrelType1,FG_COMM,req(ireq),IERR) > c write (iout,*) "ireq,req",ireq,req(ireq) > endif > enddo > C Send the contacts to processors that need them > do ii=1,ntask_cont_to > iproc=itask_cont_to(ii) > nn=ncont_sent(ii) > c write (iout,*) nn," contacts to processor",iproc, > c & " of CONT_TO_COMM group" > if (nn.gt.0) then > ireq=ireq+1 > call MPI_Isend(zapas(1,1,ii),nn*max_dim,MPI_DOUBLE_PRECISION, > & iproc,CorrelType1,FG_COMM,req(ireq),IERR) > c write (iout,*) "ireq,req",ireq,req(ireq) > c do i=1,nn > c write(iout,'(2f5.0,4f10.5)')(zapas(j,i,ii),j=1,5) > c enddo > endif > enddo > c write (iout,*) "number of requests (contacts)",ireq > c write (iout,*) "req",(req(i),i=1,4) > c call flush(iout) > if (ireq.gt.0) > & call MPI_Waitall(ireq,req,status_array,ierr) > do iii=1,ntask_cont_from > iproc=itask_cont_from(iii) > nn=ncont_recv(iii) > if (lprn) then > write (iout,*) "Received",nn," contacts from processor",iproc, > & " of CONT_FROM_COMM group" > call flush(iout) > do i=1,nn > write(iout,'(2f5.0,4f10.5)')(zapas_recv(j,i,iii),j=1,5) > enddo > call flush(iout) > endif > do i=1,nn > ii=zapas_recv(1,i,iii) > c Flag the received contacts to prevent double-counting > jj=-zapas_recv(2,i,iii) > c write (iout,*) "iii",iii," i",i," ii",ii," jj",jj > c call flush(iout) > nnn=num_cont_hb(ii)+1 > num_cont_hb(ii)=nnn > jcont_hb(nnn,ii)=jj > facont_hb(nnn,ii)=zapas_recv(3,i,iii) > ees0p(nnn,ii)=zapas_recv(4,i,iii) > ees0m(nnn,ii)=zapas_recv(5,i,iii) > gacont_hbr(1,nnn,ii)=zapas_recv(6,i,iii) > gacont_hbr(2,nnn,ii)=zapas_recv(7,i,iii) > gacont_hbr(3,nnn,ii)=zapas_recv(8,i,iii) > gacontm_hb1(1,nnn,ii)=zapas_recv(9,i,iii) > gacontm_hb1(2,nnn,ii)=zapas_recv(10,i,iii) > gacontm_hb1(3,nnn,ii)=zapas_recv(11,i,iii) > gacontp_hb1(1,nnn,ii)=zapas_recv(12,i,iii) > gacontp_hb1(2,nnn,ii)=zapas_recv(13,i,iii) > gacontp_hb1(3,nnn,ii)=zapas_recv(14,i,iii) > gacontm_hb2(1,nnn,ii)=zapas_recv(15,i,iii) > gacontm_hb2(2,nnn,ii)=zapas_recv(16,i,iii) > gacontm_hb2(3,nnn,ii)=zapas_recv(17,i,iii) > gacontp_hb2(1,nnn,ii)=zapas_recv(18,i,iii) > gacontp_hb2(2,nnn,ii)=zapas_recv(19,i,iii) > gacontp_hb2(3,nnn,ii)=zapas_recv(20,i,iii) > gacontm_hb3(1,nnn,ii)=zapas_recv(21,i,iii) > gacontm_hb3(2,nnn,ii)=zapas_recv(22,i,iii) > gacontm_hb3(3,nnn,ii)=zapas_recv(23,i,iii) > gacontp_hb3(1,nnn,ii)=zapas_recv(24,i,iii) > gacontp_hb3(2,nnn,ii)=zapas_recv(25,i,iii) > gacontp_hb3(3,nnn,ii)=zapas_recv(26,i,iii) > enddo > enddo > call flush(iout) > if (lprn) then > write (iout,'(a)') 'Contact function values after receive:' > do i=nnt,nct-2 > write (iout,'(2i3,50(1x,i3,f5.2))') > & i,num_cont_hb(i),(jcont_hb(j,i),facont_hb(j,i), > & j=1,num_cont_hb(i)) > enddo > call flush(iout) > endif 4842c6046 < write (iout,'(2i3,50(1x,i2,f5.2))') --- > write (iout,'(2i3,50(1x,i3,f5.2))') 4856c6060 < do i=iatel_s,iatel_e+1 --- > do i=min0(iatel_s,iturn4_start),max0(iatel_e,iturn3_end) 4861a6066 > jp=iabs(j) 4863a6069 > jp1=iabs(j1) 4866c6072,6074 < if (j1.eq.j+1 .or. j1.eq.j-1) then --- > if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 > & .or. j.lt.0 .and. j1.gt.0) .and. > & (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then 4869c6077,6079 < ecorr=ecorr+ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0) --- > ecorr=ecorr+ehbcorr(i,jp,i+1,jp1,jj,kk,0.72D0,0.32D0) > if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') > & 'ecorrh',i,j,ehbcorr(i,j,i+1,j1,jj,kk,0.72D0,0.32D0) 4891a6102,6158 > subroutine add_hb_contact(ii,jj,itask) > implicit real*8 (a-h,o-z) > include "DIMENSIONS" > include "COMMON.IOUNITS" > integer max_cont > integer max_dim > parameter (max_cont=maxconts) > parameter (max_dim=26) > include "COMMON.CONTACTS" > double precision zapas(max_dim,maxconts,max_fg_procs), > & zapas_recv(max_dim,maxconts,max_fg_procs) > common /przechowalnia/ zapas > integer i,j,ii,jj,iproc,itask(4),nn > c write (iout,*) "itask",itask > do i=1,2 > iproc=itask(i) > if (iproc.gt.0) then > do j=1,num_cont_hb(ii) > jjc=jcont_hb(j,ii) > c write (iout,*) "i",ii," j",jj," jjc",jjc > if (jjc.eq.jj) then > ncont_sent(iproc)=ncont_sent(iproc)+1 > nn=ncont_sent(iproc) > zapas(1,nn,iproc)=ii > zapas(2,nn,iproc)=jjc > zapas(3,nn,iproc)=facont_hb(j,ii) > zapas(4,nn,iproc)=ees0p(j,ii) > zapas(5,nn,iproc)=ees0m(j,ii) > zapas(6,nn,iproc)=gacont_hbr(1,j,ii) > zapas(7,nn,iproc)=gacont_hbr(2,j,ii) > zapas(8,nn,iproc)=gacont_hbr(3,j,ii) > zapas(9,nn,iproc)=gacontm_hb1(1,j,ii) > zapas(10,nn,iproc)=gacontm_hb1(2,j,ii) > zapas(11,nn,iproc)=gacontm_hb1(3,j,ii) > zapas(12,nn,iproc)=gacontp_hb1(1,j,ii) > zapas(13,nn,iproc)=gacontp_hb1(2,j,ii) > zapas(14,nn,iproc)=gacontp_hb1(3,j,ii) > zapas(15,nn,iproc)=gacontm_hb2(1,j,ii) > zapas(16,nn,iproc)=gacontm_hb2(2,j,ii) > zapas(17,nn,iproc)=gacontm_hb2(3,j,ii) > zapas(18,nn,iproc)=gacontp_hb2(1,j,ii) > zapas(19,nn,iproc)=gacontp_hb2(2,j,ii) > zapas(20,nn,iproc)=gacontp_hb2(3,j,ii) > zapas(21,nn,iproc)=gacontm_hb3(1,j,ii) > zapas(22,nn,iproc)=gacontm_hb3(2,j,ii) > zapas(23,nn,iproc)=gacontm_hb3(3,j,ii) > zapas(24,nn,iproc)=gacontp_hb3(1,j,ii) > zapas(25,nn,iproc)=gacontp_hb3(2,j,ii) > zapas(26,nn,iproc)=gacontp_hb3(3,j,ii) > exit > endif > enddo > endif > enddo > return > end > c------------------------------------------------------------------------------ 4897d6163 < include 'DIMENSIONS.ZSCOPT' 4899,4900c6165,6174 < #ifdef MPL < include 'COMMON.INFO' --- > #ifdef MPI > include "mpif.h" > parameter (max_cont=maxconts) > parameter (max_dim=70) > integer source,CorrelType,CorrelID,CorrelType1,CorrelID1,Error > double precision zapas(max_dim,maxconts,max_fg_procs), > & zapas_recv(max_dim,maxconts,max_fg_procs) > common /przechowalnia/ zapas > integer status(MPI_STATUS_SIZE),req(maxconts*2), > & status_array(MPI_STATUS_SIZE,maxconts*2) 4901a6176 > include 'COMMON.SETUP' 4903a6179 > include 'COMMON.LOCAL' 4906,4913c6182,6183 < #ifdef MPL < parameter (max_cont=maxconts) < parameter (max_dim=2*(8*3+2)) < parameter (msglen1=max_cont*max_dim*4) < parameter (msglen2=2*msglen1) < integer source,CorrelType,CorrelID,Error < double precision buffer(max_cont,max_dim) < #endif --- > include 'COMMON.CHAIN' > include 'COMMON.CONTROL' 4914a6185 > integer num_cont_hb_old(maxres) 4916c6187,6188 < --- > double precision eello4,eello5,eelo6,eello_turn6 > external eello4,eello5,eello6,eello_turn6 4920c6192,6195 < #ifdef MPL --- > #ifdef MPI > do i=1,nres > num_cont_hb_old(i)=num_cont_hb(i) > enddo 4923c6198 < if (fgProcs.le.1) goto 30 --- > if (nfgtasks.le.1) goto 30 4925c6200 < write (iout,'(a)') 'Contact function values:' --- > write (iout,'(a)') 'Contact function values before RECEIVE:' 4932,4933c6207,6282 < C Caution! Following code assumes that electrostatic interactions concerning < C a given atom are split among at most two processors! --- > call flush(iout) > do i=1,ntask_cont_from > ncont_recv(i)=0 > enddo > do i=1,ntask_cont_to > ncont_sent(i)=0 > enddo > c write (iout,*) "ntask_cont_from",ntask_cont_from," ntask_cont_to", > c & ntask_cont_to > C Make the list of contacts to send to send to other procesors > do i=iturn3_start,iturn3_end > c write (iout,*) "make contact list turn3",i," num_cont", > c & num_cont_hb(i) > call add_hb_contact_eello(i,i+2,iturn3_sent_local(1,i)) > enddo > do i=iturn4_start,iturn4_end > c write (iout,*) "make contact list turn4",i," num_cont", > c & num_cont_hb(i) > call add_hb_contact_eello(i,i+3,iturn4_sent_local(1,i)) > enddo > do ii=1,nat_sent > i=iat_sent(ii) > c write (iout,*) "make contact list longrange",i,ii," num_cont", > c & num_cont_hb(i) > do j=1,num_cont_hb(i) > do k=1,4 > jjc=jcont_hb(j,i) > iproc=iint_sent_local(k,jjc,ii) > c write (iout,*) "i",i," j",j," k",k," jjc",jjc," iproc",iproc > if (iproc.ne.0) then > ncont_sent(iproc)=ncont_sent(iproc)+1 > nn=ncont_sent(iproc) > zapas(1,nn,iproc)=i > zapas(2,nn,iproc)=jjc > zapas(3,nn,iproc)=d_cont(j,i) > ind=3 > do kk=1,3 > ind=ind+1 > zapas(ind,nn,iproc)=grij_hb_cont(kk,j,i) > enddo > do kk=1,2 > do ll=1,2 > ind=ind+1 > zapas(ind,nn,iproc)=a_chuj(ll,kk,j,i) > enddo > enddo > do jj=1,5 > do kk=1,3 > do ll=1,2 > do mm=1,2 > ind=ind+1 > zapas(ind,nn,iproc)=a_chuj_der(mm,ll,kk,jj,j,i) > enddo > enddo > enddo > enddo > endif > enddo > enddo > enddo > if (lprn) then > write (iout,*) > & "Numbers of contacts to be sent to other processors", > & (ncont_sent(i),i=1,ntask_cont_to) > write (iout,*) "Contacts sent" > do ii=1,ntask_cont_to > nn=ncont_sent(ii) > iproc=itask_cont_to(ii) > write (iout,*) nn," contacts to processor",iproc, > & " of CONT_TO_COMM group" > do i=1,nn > write(iout,'(2f5.0,10f10.5)')(zapas(j,i,ii),j=1,10) > enddo > enddo > call flush(iout) > endif 4935,5023c6284,6403 < CorrelID=MyID+1 < ldone=.false. < do i=1,max_cont < do j=1,max_dim < buffer(i,j)=0.0D0 < enddo < enddo < mm=mod(MyRank,2) < cd write (iout,*) 'MyRank',MyRank,' mm',mm < if (mm) 20,20,10 < 10 continue < cd write (iout,*) 'Sending: MyRank',MyRank,' mm',mm,' ldone',ldone < if (MyRank.gt.0) then < C Send correlation contributions to the preceding processor < msglen=msglen1 < nn=num_cont_hb(iatel_s) < call pack_buffer(max_cont,max_dim,iatel_s,0,buffer) < cd write (iout,*) 'The BUFFER array:' < cd do i=1,nn < cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,26) < cd enddo < if (ielstart(iatel_s).gt.iatel_s+ispp) then < msglen=msglen2 < call pack_buffer(max_cont,max_dim,iatel_s+1,26,buffer) < C Clear the contacts of the atom passed to the neighboring processor < nn=num_cont_hb(iatel_s+1) < cd do i=1,nn < cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j+26),j=1,26) < cd enddo < num_cont_hb(iatel_s)=0 < endif < cd write (iout,*) 'Processor ',MyID,MyRank, < cd & ' is sending correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen < cd write (*,*) 'Processor ',MyID,MyRank, < cd & ' is sending correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen,' CorrelType=',CorrelType < call mp_bsend(buffer,msglen,MyID-1,CorrelType,CorrelID) < cd write (iout,*) 'Processor ',MyID, < cd & ' has sent correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen,' CorrelID=',CorrelID < cd write (*,*) 'Processor ',MyID, < cd & ' has sent correlation contribution to processor',MyID-1, < cd & ' msglen=',msglen,' CorrelID=',CorrelID < msglen=msglen1 < endif ! (MyRank.gt.0) < if (ldone) goto 30 < ldone=.true. < 20 continue < cd write (iout,*) 'Receiving: MyRank',MyRank,' mm',mm,' ldone',ldone < if (MyRank.lt.fgProcs-1) then < C Receive correlation contributions from the next processor < msglen=msglen1 < if (ielend(iatel_e).lt.nct-1) msglen=msglen2 < cd write (iout,*) 'Processor',MyID, < cd & ' is receiving correlation contribution from processor',MyID+1, < cd & ' msglen=',msglen,' CorrelType=',CorrelType < cd write (*,*) 'Processor',MyID, < cd & ' is receiving correlation contribution from processor',MyID+1, < cd & ' msglen=',msglen,' CorrelType=',CorrelType < nbytes=-1 < do while (nbytes.le.0) < call mp_probe(MyID+1,CorrelType,nbytes) < enddo < cd print *,'Processor',MyID,' msglen',msglen,' nbytes',nbytes < call mp_brecv(buffer,msglen,MyID+1,CorrelType,nbytes) < cd write (iout,*) 'Processor',MyID, < cd & ' has received correlation contribution from processor',MyID+1, < cd & ' msglen=',msglen,' nbytes=',nbytes < cd write (iout,*) 'The received BUFFER array:' < cd do i=1,max_cont < cd write (iout,'(i2,9(3f8.3,2x))') i,(buffer(i,j),j=1,52) < cd enddo < if (msglen.eq.msglen1) then < call unpack_buffer(max_cont,max_dim,iatel_e+1,0,buffer) < else if (msglen.eq.msglen2) then < call unpack_buffer(max_cont,max_dim,iatel_e,0,buffer) < call unpack_buffer(max_cont,max_dim,iatel_e+1,26,buffer) < else < write (iout,*) < & 'ERROR!!!! message length changed while processing correlations.' < write (*,*) < & 'ERROR!!!! message length changed while processing correlations.' < call mp_stopall(Error) < endif ! msglen.eq.msglen1 < endif ! MyRank.lt.fgProcs-1 < if (ldone) goto 30 < ldone=.true. < goto 10 --- > CorrelID=fg_rank+1 > CorrelType1=478 > CorrelID1=nfgtasks+fg_rank+1 > ireq=0 > C Receive the numbers of needed contacts from other processors > do ii=1,ntask_cont_from > iproc=itask_cont_from(ii) > ireq=ireq+1 > call MPI_Irecv(ncont_recv(ii),1,MPI_INTEGER,iproc,CorrelType, > & FG_COMM,req(ireq),IERR) > enddo > c write (iout,*) "IRECV ended" > c call flush(iout) > C Send the number of contacts needed by other processors > do ii=1,ntask_cont_to > iproc=itask_cont_to(ii) > ireq=ireq+1 > call MPI_Isend(ncont_sent(ii),1,MPI_INTEGER,iproc,CorrelType, > & FG_COMM,req(ireq),IERR) > enddo > c write (iout,*) "ISEND ended" > c write (iout,*) "number of requests (nn)",ireq > call flush(iout) > if (ireq.gt.0) > & call MPI_Waitall(ireq,req,status_array,ierr) > c write (iout,*) > c & "Numbers of contacts to be received from other processors", > c & (ncont_recv(i),i=1,ntask_cont_from) > c call flush(iout) > C Receive contacts > ireq=0 > do ii=1,ntask_cont_from > iproc=itask_cont_from(ii) > nn=ncont_recv(ii) > c write (iout,*) "Receiving",nn," contacts from processor",iproc, > c & " of CONT_TO_COMM group" > call flush(iout) > if (nn.gt.0) then > ireq=ireq+1 > call MPI_Irecv(zapas_recv(1,1,ii),nn*max_dim, > & MPI_DOUBLE_PRECISION,iproc,CorrelType1,FG_COMM,req(ireq),IERR) > c write (iout,*) "ireq,req",ireq,req(ireq) > endif > enddo > C Send the contacts to processors that need them > do ii=1,ntask_cont_to > iproc=itask_cont_to(ii) > nn=ncont_sent(ii) > c write (iout,*) nn," contacts to processor",iproc, > c & " of CONT_TO_COMM group" > if (nn.gt.0) then > ireq=ireq+1 > call MPI_Isend(zapas(1,1,ii),nn*max_dim,MPI_DOUBLE_PRECISION, > & iproc,CorrelType1,FG_COMM,req(ireq),IERR) > c write (iout,*) "ireq,req",ireq,req(ireq) > c do i=1,nn > c write(iout,'(2f5.0,4f10.5)')(zapas(j,i,ii),j=1,5) > c enddo > endif > enddo > c write (iout,*) "number of requests (contacts)",ireq > c write (iout,*) "req",(req(i),i=1,4) > c call flush(iout) > if (ireq.gt.0) > & call MPI_Waitall(ireq,req,status_array,ierr) > do iii=1,ntask_cont_from > iproc=itask_cont_from(iii) > nn=ncont_recv(iii) > if (lprn) then > write (iout,*) "Received",nn," contacts from processor",iproc, > & " of CONT_FROM_COMM group" > call flush(iout) > do i=1,nn > write(iout,'(2f5.0,10f10.5)')(zapas_recv(j,i,iii),j=1,10) > enddo > call flush(iout) > endif > do i=1,nn > ii=zapas_recv(1,i,iii) > c Flag the received contacts to prevent double-counting > jj=-zapas_recv(2,i,iii) > c write (iout,*) "iii",iii," i",i," ii",ii," jj",jj > c call flush(iout) > nnn=num_cont_hb(ii)+1 > num_cont_hb(ii)=nnn > jcont_hb(nnn,ii)=jj > d_cont(nnn,ii)=zapas_recv(3,i,iii) > ind=3 > do kk=1,3 > ind=ind+1 > grij_hb_cont(kk,nnn,ii)=zapas_recv(ind,i,iii) > enddo > do kk=1,2 > do ll=1,2 > ind=ind+1 > a_chuj(ll,kk,nnn,ii)=zapas_recv(ind,i,iii) > enddo > enddo > do jj=1,5 > do kk=1,3 > do ll=1,2 > do mm=1,2 > ind=ind+1 > a_chuj_der(mm,ll,kk,jj,nnn,ii)=zapas_recv(ind,i,iii) > enddo > enddo > enddo > enddo > enddo > enddo > call flush(iout) > if (lprn) then > write (iout,'(a)') 'Contact function values after receive:' > do i=nnt,nct-2 > write (iout,'(2i3,50(1x,i3,5f6.3))') > & i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i), > & ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i)) > enddo > call flush(iout) > endif 5029,5031c6409,6411 < 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)) --- > write (iout,'(2i3,50(1x,i2,5f6.3))') > & i,num_cont_hb(i),(jcont_hb(j,i),d_cont(j,i), > & ((a_chuj(ll,kk,j,i),ll=1,2),kk=1,2),j=1,num_cont_hb(i)) 5049a6430 > #ifdef MOMENT 5050a6432 > #endif 5055c6437,6443 < do i=iatel_s,iatel_e+1 --- > c write (iout,*) "gradcorr5 in eello5 before loop" > c do iii=1,nres > c write (iout,'(i5,3f10.5)') > c & iii,(gradcorr5(jjj,iii),jjj=1,3) > c enddo > do i=min0(iatel_s,iturn4_start),max0(iatel_e+1,iturn3_end+1) > c write (iout,*) "corr loop i",i 5060a6449 > jp=iabs(j) 5063c6452,6453 < c write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, --- > jp1=iabs(j1) > c write (iout,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, 5065c6455,6458 < if (j1.eq.j+1 .or. j1.eq.j-1) then --- > c if (j1.eq.j+1 .or. j1.eq.j-1) then > if ((j.gt.0 .and. j1.gt.0 .or. j.gt.0 .and. j1.lt.0 > & .or. j.lt.0 .and. j1.gt.0) .and. > & (jp1.eq.jp+1 .or. jp1.eq.jp-1)) then 5075,5076c6468,6469 < c write (*,*) 'i=',i,' j=',j,' i1=',i1,' j1=',j1, < c & ' jj=',jj,' kk=',kk --- > cd write (iout,*) 'i=',i,' j=',jp,' i1=',i1,' j1=',jp1, > cd & ' jj=',jj,' kk=',kk 5085,5086c6478,6483 < cd & ' ekont=',ekont,' fprim=',fprimcont < call calc_eello(i,j,i+1,j1,jj,kk) --- > cd & ' ekont=',ekont,' fprim=',fprimcont, > cd & ' fac_prim1',fac_prim1,' fac_prim2',fac_prim2 > cd write (iout,*) "g_contij",g_contij > cd write (iout,*) "grij_hb_cont i",grij_hb_cont(:,jj,i) > cd write (iout,*) "grij_hb_cont i1",grij_hb_cont(:,jj,i1) > call calc_eello(i,jp,i+1,jp1,jj,kk) 5088c6485,6493 < & ecorr=ecorr+eello4(i,j,i+1,j1,jj,kk) --- > & ecorr=ecorr+eello4(i,jp,i+1,jp1,jj,kk) > if (energy_dec.and.wcorr4.gt.0.0d0) > 1 write (iout,'(a6,4i5,0pf7.3)') > 2 'ecorr4',i,j,i+1,j1,eello4(i,jp,i+1,jp1,jj,kk) > c write (iout,*) "gradcorr5 before eello5" > c do iii=1,nres > c write (iout,'(i5,3f10.5)') > c & iii,(gradcorr5(jjj,iii),jjj=1,3) > c enddo 5090,5091c6495,6503 < & ecorr5=ecorr5+eello5(i,j,i+1,j1,jj,kk) < c print *,"wcorr5",ecorr5 --- > & ecorr5=ecorr5+eello5(i,jp,i+1,jp1,jj,kk) > c write (iout,*) "gradcorr5 after eello5" > c do iii=1,nres > c write (iout,'(i5,3f10.5)') > c & iii,(gradcorr5(jjj,iii),jjj=1,3) > c enddo > if (energy_dec.and.wcorr5.gt.0.0d0) > 1 write (iout,'(a6,4i5,0pf7.3)') > 2 'ecorr5',i,j,i+1,j1,eello5(i,jp,i+1,jp1,jj,kk) 5093,5094c6505,6506 < cd write(2,*)'ijkl',i,j,i+1,j1 < if (wcorr6.gt.0.0d0 .and. (j.ne.i+4 .or. j1.ne.i+3 --- > cd write(2,*)'ijkl',i,jp,i+1,jp1 > if (wcorr6.gt.0.0d0 .and. (jp.ne.i+4 .or. jp1.ne.i+3 5097c6509,6511 < ecorr6=ecorr6+eello6(i,j,i+1,j1,jj,kk) --- > ecorr6=ecorr6+eello6(i,jp,i+1,jp1,jj,kk) > if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') > 1 'ecorr6',i,j,i+1,j1,eello6(i,jp,i+1,jp1,jj,kk) 5101,5103c6515,6517 < cd & dabs(eello4(i,j,i+1,j1,jj,kk)), < cd & dabs(eello5(i,j,i+1,j1,jj,kk)), < cd & dabs(eello6(i,j,i+1,j1,jj,kk)) --- > cd & dabs(eello4(i,jp,i+1,jp1,jj,kk)), > cd & dabs(eello5(i,jp,i+1,jp1,jj,kk)), > cd & dabs(eello6(i,jp,i+1,jp1,jj,kk)) 5105,5106c6519,6520 < & .and. (j.eq.i+4 .and. j1.eq.i+3)) then < cd write (iout,*) '******eturn6: i,j,i+1,j1',i,j,i+1,j1 --- > & .and. (jp.eq.i+4 .and. jp1.eq.i+3)) then > cd write (iout,*) '******eturn6: i,j,i+1,j1',i,jip,i+1,jp1 5107a6522,6523 > if (energy_dec) write (iout,'(a6,4i5,0pf7.3)') > 1 'eturn6',i,j,i+1,j1,eello_turn6(i,jj,kk) 5112,5115d6527 < else if (j1.eq.j) then < C Contacts I-J and I-(J+1) occur simultaneously. < C The system loses extra energy. < c ecorr=ecorr+ehbcorr(i,j,i+1,j,jj,kk,0.60D0,-0.40D0) 5118,5127d6529 < do kk=1,num_conti < j1=jcont_hb(kk,i) < c write (iout,*) '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. < c ecorr=ecorr+ehbcorr(i,j,i,j+1,jj,kk,0.60D0,-0.40D0) < endif ! j1==j+1 < enddo ! kk 5129a6532,6594 > do i=1,nres > num_cont_hb(i)=num_cont_hb_old(i) > enddo > c write (iout,*) "gradcorr5 in eello5" > c do iii=1,nres > c write (iout,'(i5,3f10.5)') > c & iii,(gradcorr5(jjj,iii),jjj=1,3) > c enddo > return > end > c------------------------------------------------------------------------------ > subroutine add_hb_contact_eello(ii,jj,itask) > implicit real*8 (a-h,o-z) > include "DIMENSIONS" > include "COMMON.IOUNITS" > integer max_cont > integer max_dim > parameter (max_cont=maxconts) > parameter (max_dim=70) > include "COMMON.CONTACTS" > double precision zapas(max_dim,maxconts,max_fg_procs), > & zapas_recv(max_dim,maxconts,max_fg_procs) > common /przechowalnia/ zapas > integer i,j,ii,jj,iproc,itask(4),nn > c write (iout,*) "itask",itask > do i=1,2 > iproc=itask(i) > if (iproc.gt.0) then > do j=1,num_cont_hb(ii) > jjc=jcont_hb(j,ii) > c write (iout,*) "send turns i",ii," j",jj," jjc",jjc > if (jjc.eq.jj) then > ncont_sent(iproc)=ncont_sent(iproc)+1 > nn=ncont_sent(iproc) > zapas(1,nn,iproc)=ii > zapas(2,nn,iproc)=jjc > zapas(3,nn,iproc)=d_cont(j,ii) > ind=3 > do kk=1,3 > ind=ind+1 > zapas(ind,nn,iproc)=grij_hb_cont(kk,j,ii) > enddo > do kk=1,2 > do ll=1,2 > ind=ind+1 > zapas(ind,nn,iproc)=a_chuj(ll,kk,j,ii) > enddo > enddo > do jj=1,5 > do kk=1,3 > do ll=1,2 > do mm=1,2 > ind=ind+1 > zapas(ind,nn,iproc)=a_chuj_der(mm,ll,kk,jj,j,ii) > enddo > enddo > enddo > enddo > exit > endif > enddo > endif > enddo 5157,5161c6622,6626 < c write (iout,*)'Contacts have occurred for peptide groups',i,j, < c & ' and',k,l < c write (iout,*)'Contacts have occurred for peptide groups', < c & i,j,' fcont:',eij,' eij',' eesij',ees0pij,ees0mij,' and ',k,l < c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' ees=',ees --- > c write (iout,'(2(a,2i3,a,f10.5,a,2f10.5),a,f10.5,a,$)') > c & 'Contacts ',i,j, > c & ' eij',eij,' eesij',ees0pij,ees0mij,' and ',k,l > c & ,' fcont ',ekl,' eeskl',ees0pkl,ees0mkl,' energy=',ekont*ees, > c & 'gradcorr_long' 5163,5164c6628 < ecorr=ecorr+ekont*ees < if (calc_grad) then --- > c ecorr=ecorr+ekont*ees 5165a6630,6633 > coeffpees0pij=coeffp*ees0pij > coeffmees0mij=coeffm*ees0mij > coeffpees0pkl=coeffp*ees0pkl > coeffmees0mkl=coeffm*ees0mkl 5167,5198c6635,6678 < ghalf=0.5D0*ees*ekl*gacont_hbr(ll,jj,i) < gradcorr(ll,i)=gradcorr(ll,i)+ghalf < & -ekont*(coeffp*ees0pkl*gacontp_hb1(ll,jj,i)+ < & coeffm*ees0mkl*gacontm_hb1(ll,jj,i)) < gradcorr(ll,j)=gradcorr(ll,j)+ghalf < & -ekont*(coeffp*ees0pkl*gacontp_hb2(ll,jj,i)+ < & coeffm*ees0mkl*gacontm_hb2(ll,jj,i)) < ghalf=0.5D0*ees*eij*gacont_hbr(ll,kk,k) < gradcorr(ll,k)=gradcorr(ll,k)+ghalf < & -ekont*(coeffp*ees0pij*gacontp_hb1(ll,kk,k)+ < & coeffm*ees0mij*gacontm_hb1(ll,kk,k)) < gradcorr(ll,l)=gradcorr(ll,l)+ghalf < & -ekont*(coeffp*ees0pij*gacontp_hb2(ll,kk,k)+ < & coeffm*ees0mij*gacontm_hb2(ll,kk,k)) < enddo < do m=i+1,j-1 < do ll=1,3 < gradcorr(ll,m)=gradcorr(ll,m)+ < & ees*ekl*gacont_hbr(ll,jj,i)- < & ekont*(coeffp*ees0pkl*gacontp_hb3(ll,jj,i)+ < & coeffm*ees0mkl*gacontm_hb3(ll,jj,i)) < enddo < enddo < do m=k+1,l-1 < do ll=1,3 < gradcorr(ll,m)=gradcorr(ll,m)+ < & ees*eij*gacont_hbr(ll,kk,k)- < & ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+ < & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) < enddo < enddo < endif --- > cgrad ghalfi=ees*ekl*gacont_hbr(ll,jj,i) > gradcorr(ll,i)=gradcorr(ll,i)!+0.5d0*ghalfi > & -ekont*(coeffpees0pkl*gacontp_hb1(ll,jj,i)+ > & coeffmees0mkl*gacontm_hb1(ll,jj,i)) > gradcorr(ll,j)=gradcorr(ll,j)!+0.5d0*ghalfi > & -ekont*(coeffpees0pkl*gacontp_hb2(ll,jj,i)+ > & coeffmees0mkl*gacontm_hb2(ll,jj,i)) > cgrad ghalfk=ees*eij*gacont_hbr(ll,kk,k) > gradcorr(ll,k)=gradcorr(ll,k)!+0.5d0*ghalfk > & -ekont*(coeffpees0pij*gacontp_hb1(ll,kk,k)+ > & coeffmees0mij*gacontm_hb1(ll,kk,k)) > gradcorr(ll,l)=gradcorr(ll,l)!+0.5d0*ghalfk > & -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_long(ll,j)=gradcorr_long(ll,j)+gradlongij > gradcorr_long(ll,i)=gradcorr_long(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_long(ll,l)=gradcorr_long(ll,l)+gradlongkl > gradcorr_long(ll,k)=gradcorr_long(ll,k)-gradlongkl > c write (iout,'(2f10.5,2x,$)') gradlongij,gradlongkl > enddo > c write (iout,*) > cgrad do m=i+1,j-1 > cgrad do ll=1,3 > cgrad gradcorr(ll,m)=gradcorr(ll,m)+ > cgrad & ees*ekl*gacont_hbr(ll,jj,i)- > cgrad & ekont*(coeffp*ees0pkl*gacontp_hb3(ll,jj,i)+ > cgrad & coeffm*ees0mkl*gacontm_hb3(ll,jj,i)) > cgrad enddo > cgrad enddo > cgrad do m=k+1,l-1 > cgrad do ll=1,3 > cgrad gradcorr(ll,m)=gradcorr(ll,m)+ > cgrad & ees*eij*gacont_hbr(ll,kk,k)- > cgrad & ekont*(coeffp*ees0pij*gacontp_hb3(ll,kk,k)+ > cgrad & coeffm*ees0mij*gacontm_hb3(ll,kk,k)) > cgrad enddo > cgrad enddo > c write (iout,*) "ehbcorr",ekont*ees 5201a6682 > #ifdef MOMENT 5206d6686 < include 'DIMENSIONS.ZSCOPT' 5240d6719 < if (.not.calc_grad) return 5264a6744 > #endif 5273d6752 < include 'DIMENSIONS.ZSCOPT' 5289a6769,6770 > cd write (iout,*) "a_chujij",((a_chuj(iii,jjj,jj,i),iii=1,2),jjj=1,2) > cd write (iout,*) "a_chujkl",((a_chuj(iii,jjj,kk,k),iii=1,2),jjj=1,2) 5650d7130 < include 'DIMENSIONS.ZSCOPT' 5671d7150 < if (calc_grad) then 5713,5718c7192,7197 < cold ghalf=0.5d0*eel4*ekl*gacont_hbr(ll,jj,i) < ggg1(ll)=eel4*g_contij(ll,1) < ggg2(ll)=eel4*g_contij(ll,2) < ghalf=0.5d0*ggg1(ll) < cd ghalf=0.0d0 < gradcorr(ll,i)=gradcorr(ll,i)+ghalf+ekont*derx(ll,2,1) --- > cgrad ggg1(ll)=eel4*g_contij(ll,1) > cgrad ggg2(ll)=eel4*g_contij(ll,2) > glongij=eel4*g_contij(ll,1)+ekont*derx(ll,1,1) > glongkl=eel4*g_contij(ll,2)+ekont*derx(ll,1,2) > cgrad ghalf=0.5d0*ggg1(ll) > gradcorr(ll,i)=gradcorr(ll,i)+ekont*derx(ll,2,1) 5720c7199 < gradcorr(ll,j)=gradcorr(ll,j)+ghalf+ekont*derx(ll,4,1) --- > gradcorr(ll,j)=gradcorr(ll,j)+ekont*derx(ll,4,1) 5722,5725c7201,7204 < cold ghalf=0.5d0*eel4*eij*gacont_hbr(ll,kk,k) < ghalf=0.5d0*ggg2(ll) < cd ghalf=0.0d0 < gradcorr(ll,k)=gradcorr(ll,k)+ghalf+ekont*derx(ll,2,2) --- > gradcorr_long(ll,j)=gradcorr_long(ll,j)+glongij > gradcorr_long(ll,i)=gradcorr_long(ll,i)-glongij > cgrad ghalf=0.5d0*ggg2(ll) > gradcorr(ll,k)=gradcorr(ll,k)+ekont*derx(ll,2,2) 5727c7206 < gradcorr(ll,l)=gradcorr(ll,l)+ghalf+ekont*derx(ll,4,2) --- > gradcorr(ll,l)=gradcorr(ll,l)+ekont*derx(ll,4,2) 5728a7208,7209 > gradcorr_long(ll,l)=gradcorr_long(ll,l)+glongkl > gradcorr_long(ll,k)=gradcorr_long(ll,k)-glongkl 5730,5753c7211,7230 < cd goto 1112 < do m=i+1,j-1 < do ll=1,3 < cold gradcorr(ll,m)=gradcorr(ll,m)+eel4*ekl*gacont_hbr(ll,jj,i) < gradcorr(ll,m)=gradcorr(ll,m)+ggg1(ll) < enddo < enddo < do m=k+1,l-1 < do ll=1,3 < cold gradcorr(ll,m)=gradcorr(ll,m)+eel4*eij*gacont_hbr(ll,kk,k) < gradcorr(ll,m)=gradcorr(ll,m)+ggg2(ll) < enddo < enddo < 1112 continue < do m=i+2,j2 < do ll=1,3 < gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,1) < enddo < enddo < do m=k+2,l2 < do ll=1,3 < gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,2) < enddo < enddo --- > cgrad do m=i+1,j-1 > cgrad do ll=1,3 > cgrad gradcorr(ll,m)=gradcorr(ll,m)+ggg1(ll) > cgrad enddo > cgrad enddo > cgrad do m=k+1,l-1 > cgrad do ll=1,3 > cgrad gradcorr(ll,m)=gradcorr(ll,m)+ggg2(ll) > cgrad enddo > cgrad enddo > cgrad do m=i+2,j2 > cgrad do ll=1,3 > cgrad gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,1) > cgrad enddo > cgrad enddo > cgrad do m=k+2,l2 > cgrad do ll=1,3 > cgrad gradcorr(ll,m)=gradcorr(ll,m)+ekont*derx(ll,1,2) > cgrad enddo > cgrad enddo 5757d7233 < endif 5767d7242 < include 'DIMENSIONS.ZSCOPT' 5847d7321 < if (calc_grad) then 5886d7359 < endif 5895d7367 < if (calc_grad) then 5926d7397 < endif 5938d7408 < if (calc_grad) then 5971d7440 < endif 5980d7448 < if (calc_grad) then 6004d7471 < endif 6015d7481 < if (calc_grad) then 6048d7513 < endif 6057d7521 < if (calc_grad) then 6082d7545 < endif 6094d7556 < if (calc_grad) then 6112a7575,7578 > C 2/11/08 AL Gradients over DC's connecting interacting sites will be > C summed up outside the subrouine as for the other subroutines > C handling long-range interactions. The old code is commented out > C with "cgrad" to keep track of changes. 6114,6115c7580,7591 < ggg1(ll)=eel5*g_contij(ll,1) < ggg2(ll)=eel5*g_contij(ll,2) --- > cgrad ggg1(ll)=eel5*g_contij(ll,1) > cgrad ggg2(ll)=eel5*g_contij(ll,2) > gradcorr5ij=eel5*g_contij(ll,1)+ekont*derx(ll,1,1) > gradcorr5kl=eel5*g_contij(ll,2)+ekont*derx(ll,1,2) > c write (iout,'(a,3i3,a,5f8.3,2i3,a,5f8.3,a,f8.3)') > c & "ecorr5",ll,i,j," derx",derx(ll,2,1),derx(ll,3,1),derx(ll,4,1), > c & derx(ll,5,1),k,l," derx",derx(ll,2,2),derx(ll,3,2), > c & derx(ll,4,2),derx(ll,5,2)," ekont",ekont > c write (iout,'(a,3i3,a,3f8.3,2i3,a,3f8.3)') > c & "ecorr5",ll,i,j," gradcorr5",g_contij(ll,1),derx(ll,1,1), > c & gradcorr5ij, > c & k,l," gradcorr5",g_contij(ll,2),derx(ll,1,2),gradcorr5kl 6117c7593 < ghalf=0.5d0*ggg1(ll) --- > cgrad ghalf=0.5d0*ggg1(ll) 6119c7595 < gradcorr5(ll,i)=gradcorr5(ll,i)+ghalf+ekont*derx(ll,2,1) --- > gradcorr5(ll,i)=gradcorr5(ll,i)+ekont*derx(ll,2,1) 6121c7597 < gradcorr5(ll,j)=gradcorr5(ll,j)+ghalf+ekont*derx(ll,4,1) --- > gradcorr5(ll,j)=gradcorr5(ll,j)+ekont*derx(ll,4,1) 6122a7599,7600 > gradcorr5_long(ll,j)=gradcorr5_long(ll,j)+gradcorr5ij > gradcorr5_long(ll,i)=gradcorr5_long(ll,i)-gradcorr5ij 6124c7602 < ghalf=0.5d0*ggg2(ll) --- > cgrad ghalf=0.5d0*ggg2(ll) 6129a7608,7609 > gradcorr5_long(ll,l)=gradcorr5_long(ll,l)+gradcorr5kl > gradcorr5_long(ll,k)=gradcorr5_long(ll,k)-gradcorr5kl 6132,6133c7612,7613 < do m=i+1,j-1 < do ll=1,3 --- > cgrad do m=i+1,j-1 > cgrad do ll=1,3 6135,6139c7615,7619 < gradcorr5(ll,m)=gradcorr5(ll,m)+ggg1(ll) < enddo < enddo < do m=k+1,l-1 < do ll=1,3 --- > cgrad gradcorr5(ll,m)=gradcorr5(ll,m)+ggg1(ll) > cgrad enddo > cgrad enddo > cgrad do m=k+1,l-1 > cgrad do ll=1,3 6141,6143c7621,7623 < gradcorr5(ll,m)=gradcorr5(ll,m)+ggg2(ll) < enddo < enddo --- > cgrad gradcorr5(ll,m)=gradcorr5(ll,m)+ggg2(ll) > cgrad enddo > cgrad enddo 6145,6154c7625,7634 < do m=i+2,j2 < do ll=1,3 < gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,1) < enddo < enddo < do m=k+2,l2 < do ll=1,3 < gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,2) < enddo < enddo --- > cgrad do m=i+2,j2 > cgrad do ll=1,3 > cgrad gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,1) > cgrad enddo > cgrad enddo > cgrad do m=k+2,l2 > cgrad do ll=1,3 > cgrad gradcorr5(ll,m)=gradcorr5(ll,m)+ekont*derx(ll,1,2) > cgrad enddo > cgrad enddo 6158d7637 < endif 6168d7646 < include 'DIMENSIONS.ZSCOPT' 6228,6233c7706,7711 < cd write(iout,*) 'eello6_1',eello6_1,' eel6_1_num',16*eel6_1_num < cd write(iout,*) 'eello6_2',eello6_2,' eel6_2_num',16*eel6_2_num < cd write(iout,*) 'eello6_3',eello6_3,' eel6_3_num',16*eel6_3_num < cd write(iout,*) 'eello6_4',eello6_4,' eel6_4_num',16*eel6_4_num < cd write(iout,*) 'eello6_5',eello6_5,' eel6_5_num',16*eel6_5_num < cd write(iout,*) 'eello6_6',eello6_6,' eel6_6_num',16*eel6_6_num --- > cd write(iout,*) 'eello6_1',eello6_1!,' eel6_1_num',16*eel6_1_num > cd write(iout,*) 'eello6_2',eello6_2!,' eel6_2_num',16*eel6_2_num > cd write(iout,*) 'eello6_3',eello6_3!,' eel6_3_num',16*eel6_3_num > cd write(iout,*) 'eello6_4',eello6_4!,' eel6_4_num',16*eel6_4_num > cd write(iout,*) 'eello6_5',eello6_5!,' eel6_5_num',16*eel6_5_num > cd write(iout,*) 'eello6_6',eello6_6!,' eel6_6_num',16*eel6_6_num 6235d7712 < if (calc_grad) then 6251,6252c7728,7729 < ggg1(ll)=eel6*g_contij(ll,1) < ggg2(ll)=eel6*g_contij(ll,2) --- > cgrad ggg1(ll)=eel6*g_contij(ll,1) > cgrad ggg2(ll)=eel6*g_contij(ll,2) 6254c7731 < ghalf=0.5d0*ggg1(ll) --- > cgrad ghalf=0.5d0*ggg1(ll) 6256c7733,7735 < gradcorr6(ll,i)=gradcorr6(ll,i)+ghalf+ekont*derx(ll,2,1) --- > gradcorr6ij=eel6*g_contij(ll,1)+ekont*derx(ll,1,1) > gradcorr6kl=eel6*g_contij(ll,2)+ekont*derx(ll,1,2) > gradcorr6(ll,i)=gradcorr6(ll,i)+ekont*derx(ll,2,1) 6258c7737 < gradcorr6(ll,j)=gradcorr6(ll,j)+ghalf+ekont*derx(ll,4,1) --- > gradcorr6(ll,j)=gradcorr6(ll,j)+ekont*derx(ll,4,1) 6260c7739,7741 < ghalf=0.5d0*ggg2(ll) --- > gradcorr6_long(ll,j)=gradcorr6_long(ll,j)+gradcorr6ij > gradcorr6_long(ll,i)=gradcorr6_long(ll,i)-gradcorr6ij > cgrad ghalf=0.5d0*ggg2(ll) 6263c7744 < gradcorr6(ll,k)=gradcorr6(ll,k)+ghalf+ekont*derx(ll,2,2) --- > gradcorr6(ll,k)=gradcorr6(ll,k)+ekont*derx(ll,2,2) 6265c7746 < gradcorr6(ll,l)=gradcorr6(ll,l)+ghalf+ekont*derx(ll,4,2) --- > gradcorr6(ll,l)=gradcorr6(ll,l)+ekont*derx(ll,4,2) 6266a7748,7749 > gradcorr6_long(ll,l)=gradcorr6_long(ll,l)+gradcorr6kl > gradcorr6_long(ll,k)=gradcorr6_long(ll,k)-gradcorr6kl 6269,6270c7752,7753 < do m=i+1,j-1 < do ll=1,3 --- > cgrad do m=i+1,j-1 > cgrad do ll=1,3 6272,6276c7755,7759 < gradcorr6(ll,m)=gradcorr6(ll,m)+ggg1(ll) < enddo < enddo < do m=k+1,l-1 < do ll=1,3 --- > cgrad gradcorr6(ll,m)=gradcorr6(ll,m)+ggg1(ll) > cgrad enddo > cgrad enddo > cgrad do m=k+1,l-1 > cgrad do ll=1,3 6278,6291c7761,7774 < gradcorr6(ll,m)=gradcorr6(ll,m)+ggg2(ll) < enddo < enddo < 1112 continue < do m=i+2,j2 < do ll=1,3 < gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,1) < enddo < enddo < do m=k+2,l2 < do ll=1,3 < gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,2) < enddo < enddo --- > cgrad gradcorr6(ll,m)=gradcorr6(ll,m)+ggg2(ll) > cgrad enddo > cgrad enddo > cgrad1112 continue > cgrad do m=i+2,j2 > cgrad do ll=1,3 > cgrad gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,1) > cgrad enddo > cgrad enddo > cgrad do m=k+2,l2 > cgrad do ll=1,3 > cgrad gradcorr6(ll,m)=gradcorr6(ll,m)+ekont*derx(ll,1,2) > cgrad enddo > cgrad enddo 6295d7777 < endif 6305d7786 < include 'DIMENSIONS.ZSCOPT' 6346d7826 < if (.not. calc_grad) return 6411d7890 < include 'DIMENSIONS.ZSCOPT' 6461d7939 < if (.not. calc_grad) return 6598d8075 < include 'DIMENSIONS.ZSCOPT' 6651c8128,8129 < cd write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4 --- > cd write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4, > cd & "sum",-(s2+s3+s4) 6658d8135 < if (.not. calc_grad) return 6714d8190 < include 'DIMENSIONS.ZSCOPT' 6794d8269 < if (.not. calc_grad) return 6960d8434 < include 'DIMENSIONS.ZSCOPT' 6975a8450,8453 > s1=0.0d0 > s8=0.0d0 > s13=0.0d0 > c 7013,7014d8490 < #else < s1 = 0.0d0 7024,7025d8499 < #else < s8=0.0d0 7037,7038d8510 < #else < s13=0.0d0 7047d8518 < if (calc_grad) then 7048a8520,8521 > s1d =0.0d0 > s8d =0.0d0 7056,7057d8528 < #else < s8d=0.0d0 7074,7075d8544 < #else < s1d=0.0d0 7089,7090d8557 < #else < s13d=0.0d0 7112,7113d8578 < #else < s13d = 0.0d0 7130,7131d8594 < #else < s1d = 0.0d0 7140,7141d8602 < #else < s8d = 0.0d0 7149,7150d8609 < #else < s13d = 0.0d0 7172,7173d8630 < #else < s1d = 0.0d0 7184,7185d8640 < #else < s8d = 0.0d0 7248,7250c8703,8705 < ggg1(ll)=eel_turn6*g_contij(ll,1) < ggg2(ll)=eel_turn6*g_contij(ll,2) < ghalf=0.5d0*ggg1(ll) --- > cgrad ggg1(ll)=eel_turn6*g_contij(ll,1) > cgrad ggg2(ll)=eel_turn6*g_contij(ll,2) > cgrad ghalf=0.5d0*ggg1(ll) 7252c8707,8709 < gcorr6_turn(ll,i)=gcorr6_turn(ll,i)+ghalf --- > gturn6ij=eel_turn6*g_contij(ll,1)+ekont*derx_turn(ll,1,1) > gturn6kl=eel_turn6*g_contij(ll,2)+ekont*derx_turn(ll,1,2) > gcorr6_turn(ll,i)=gcorr6_turn(ll,i)!+ghalf 7255c8712 < gcorr6_turn(ll,j)=gcorr6_turn(ll,j)+ghalf --- > gcorr6_turn(ll,j)=gcorr6_turn(ll,j)!+ghalf 7258c8715,8717 < ghalf=0.5d0*ggg2(ll) --- > gcorr6_turn_long(ll,j)=gcorr6_turn_long(ll,j)+gturn6ij > gcorr6_turn_long(ll,i)=gcorr6_turn_long(ll,i)-gturn6ij > cgrad ghalf=0.5d0*ggg2(ll) 7260c8719 < gcorr6_turn(ll,k)=gcorr6_turn(ll,k)+ghalf --- > gcorr6_turn(ll,k)=gcorr6_turn(ll,k)!+ghalf 7263c8722 < gcorr6_turn(ll,l)=gcorr6_turn(ll,l)+ghalf --- > gcorr6_turn(ll,l)=gcorr6_turn(ll,l)!+ghalf 7265a8725,8726 > gcorr6_turn_long(ll,l)=gcorr6_turn_long(ll,l)+gturn6kl > gcorr6_turn_long(ll,k)=gcorr6_turn_long(ll,k)-gturn6kl 7268,7288c8729,8749 < do m=i+1,j-1 < do ll=1,3 < gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg1(ll) < enddo < enddo < do m=k+1,l-1 < do ll=1,3 < gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg2(ll) < enddo < enddo < 1112 continue < do m=i+2,j2 < do ll=1,3 < gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,1) < enddo < enddo < do m=k+2,l2 < do ll=1,3 < gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,2) < enddo < enddo --- > cgrad do m=i+1,j-1 > cgrad do ll=1,3 > cgrad gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg1(ll) > cgrad enddo > cgrad enddo > cgrad do m=k+1,l-1 > cgrad do ll=1,3 > cgrad gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ggg2(ll) > cgrad enddo > cgrad enddo > cgrad1112 continue > cgrad do m=i+2,j2 > cgrad do ll=1,3 > cgrad gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,1) > cgrad enddo > cgrad enddo > cgrad do m=k+2,l2 > cgrad do ll=1,3 > cgrad gcorr6_turn(ll,m)=gcorr6_turn(ll,m)+ekont*derx_turn(ll,1,2) > cgrad enddo > cgrad enddo 7292d8752 < endif 7297a8758,8777 > > C----------------------------------------------------------------------------- > double precision function scalar(u,v) > !DIR$ INLINEALWAYS scalar > #ifndef OSF > cDEC$ ATTRIBUTES FORCEINLINE::scalar > #endif > implicit none > double precision u(3),v(3) > cd double precision sc > cd integer i > cd sc=0.0d0 > cd do i=1,3 > cd sc=sc+u(i)*v(i) > cd enddo > cd scalar=sc > > scalar=u(1)*v(1)+u(2)*v(2)+u(3)*v(3) > return > end 7299a8780,8783 > !DIR$ INLINEALWAYS MATVEC2 > #ifndef OSF > cDEC$ ATTRIBUTES FORCEINLINE::MATVEC2 > #endif 7317a8802,8804 > #ifndef OSF > cDEC$ ATTRIBUTES FORCEINLINE::MATMAT2 > #endif 7343a8831 > !DIR$ INLINEALWAYS scalar2 7354a8843,8846 > !DIR$ INLINEALWAYS transpose2 > #ifndef OSF > cDEC$ ATTRIBUTES FORCEINLINE::transpose2 > #endif 7376a8869,8872 > !DIR$ INLINEALWAYS prodmat3 > #ifndef OSF > cDEC$ ATTRIBUTES FORCEINLINE::prodmat3 > #endif 7419,7431d8914 < C----------------------------------------------------------------------------- < double precision function scalar(u,v) < implicit none < double precision u(3),v(3) < double precision sc < integer i < sc=0.0d0 < do i=1,3 < sc=sc+u(i)*v(i) < enddo < scalar=sc < return < end