+++ /dev/null
-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<i'
-> 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
projects / unres.git / blobdiff