X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fdiffer;fp=source%2Fwham%2Fsrc-M%2Fdiffer;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=48384a1cf7705ef74cd878417d203583ebf026bd;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/wham/src-M/differ b/source/wham/src-M/differ deleted file mode 100644 index 48384a1..0000000 --- a/source/wham/src-M/differ +++ /dev/null @@ -1,5173 +0,0 @@ -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