+! 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
+! write (iout,*),"Processor",myrank," BROADCAST weights"
+ call MPI_Bcast(c(1,1),nres6,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST c"
+ call MPI_Bcast(dc(1,1),nres6,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST dc"
+ call MPI_Bcast(dc_norm(1,1),nres6,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST dc_norm"
+ call MPI_Bcast(theta(1),nres,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST theta"
+ call MPI_Bcast(phi(1),nres,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST phi"
+ call MPI_Bcast(alph(1),nres,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST alph"
+ call MPI_Bcast(omeg(1),nres,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST omeg"
+ call MPI_Bcast(vbld(1),2*nres,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+! write (iout,*) "Processor",myrank," BROADCAST vbld"
+ call MPI_Bcast(vbld_inv(1),2*nres,MPI_DOUBLE_PRECISION,&
+ king,FG_COMM,IERR)
+ time_Bcast=time_Bcast+MPI_Wtime()-time00
+! write (iout,*) "Processor",myrank," BROADCAST vbld_inv"
+ endif
+! write (iout,*) 'Processor',myrank,
+! & ' calling etotal_short ipot=',ipot
+! call flush(iout)
+! print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
+#endif
+! call int_from_cart1(.false.)
+!
+! Compute the side-chain and electrostatic interaction energy
+!
+ goto (101,102,103,104,105,106) ipot
+! Lennard-Jones potential.
+ 101 call elj_short(evdw)
+!d print '(a)','Exit ELJ'
+ goto 107
+! Lennard-Jones-Kihara potential (shifted).
+ 102 call eljk_short(evdw)
+ goto 107
+! Berne-Pechukas potential (dilated LJ, angular dependence).
+ 103 call ebp_short(evdw)
+ goto 107
+! Gay-Berne potential (shifted LJ, angular dependence).
+ 104 call egb_short(evdw)
+ goto 107
+! Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
+ 105 call egbv_short(evdw)
+ goto 107
+! Soft-sphere potential - already dealt with in the long-range part
+ 106 evdw=0.0d0
+! 106 call e_softsphere_short(evdw)
+!
+! Calculate electrostatic (H-bonding) energy of the main chain.
+!
+ 107 continue
+!
+! Calculate the short-range part of Evdwpp
+!
+ call evdwpp_short(evdw1)
+!
+! Calculate the short-range part of ESCp
+!
+ if (ipot.lt.6) then
+ call escp_short(evdw2,evdw2_14)
+ endif
+!
+! Calculate the bond-stretching energy
+!
+ call ebond(estr)
+!
+! Calculate the disulfide-bridge and other energy and the contributions
+! from other distance constraints.
+! call edis(ehpb)
+!
+! Calculate the virtual-bond-angle energy.
+!
+! Calculate the SC local energy.
+!
+ call vec_and_deriv
+ call esc(escloc)
+!
+ if (wang.gt.0d0) then
+ if (tor_mode.eq.0) then
+ call ebend(ebe)
+ else
+!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+ call ebend_kcc(ebe)
+ endif
+ else
+ ebe=0.0d0
+ endif
+ ethetacnstr=0.0d0
+ if (with_theta_constr) call etheta_constr(ethetacnstr)
+
+! write(iout,*) "in etotal afer ebe",ipot
+
+! print *,"Processor",myrank," computed UB"
+!
+! Calculate the SC local energy.
+!
+ call esc(escloc)
+!elwrite(iout,*) "in etotal afer esc",ipot
+! print *,"Processor",myrank," computed USC"
+!
+! Calculate the virtual-bond torsional energy.
+!
+!d print *,'nterm=',nterm
+! if (wtor.gt.0) then
+! call etor(etors,edihcnstr)
+! else
+! etors=0
+! edihcnstr=0
+! endif
+ if (wtor.gt.0.0d0) then
+ if (tor_mode.eq.0) then
+ call etor(etors)
+ else
+!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the
+!C energy function
+ call etor_kcc(etors)
+ endif
+ else
+ etors=0.0d0
+ endif
+ edihcnstr=0.0d0
+ if (ndih_constr.gt.0) call etor_constr(edihcnstr)
+
+! Calculate the virtual-bond torsional energy.
+!
+!
+! 6/23/01 Calculate double-torsional energy
+!
+ if ((wtor_d.gt.0.0d0).and.(tor_mode.eq.0)) then
+ call etor_d(etors_d)
+ endif
+!
+! Homology restraints
+!
+ if (constr_homology.ge.1) then
+ call e_modeller(ehomology_constr)
+! print *,"tu"
+ else
+ ehomology_constr=0.0d0
+ endif
+
+!
+! 21/5/07 Calculate local sicdechain correlation energy
+!
+ if (wsccor.gt.0.0d0) then
+ call eback_sc_corr(esccor)
+ else
+ esccor=0.0d0
+ endif
+!
+! Put energy components into an array
+!
+ do i=1,n_ene
+ energia(i)=0.0d0
+ enddo
+ energia(1)=evdw
+#ifdef SCP14
+ energia(2)=evdw2-evdw2_14
+ energia(18)=evdw2_14
+#else
+ energia(2)=evdw2
+ energia(18)=0.0d0
+#endif
+#ifdef SPLITELE
+ energia(16)=evdw1
+#else
+ energia(3)=evdw1
+#endif
+ energia(11)=ebe
+ energia(12)=escloc
+ energia(13)=etors
+ energia(14)=etors_d
+ energia(15)=ehpb
+ energia(17)=estr
+ energia(19)=edihcnstr
+ energia(21)=esccor
+ energia(51)=ehomology_constr
+! write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
+ call flush(iout)
+ call sum_energy(energia,.true.)
+! write (iout,*) "Exit ETOTAL_SHORT"
+ call flush(iout)
+ return
+ end subroutine etotal_short
+!-----------------------------------------------------------------------------
+! gnmr1.f
+!-----------------------------------------------------------------------------
+ real(kind=8) function gnmr1(y,ymin,ymax)
+! implicit none
+ real(kind=8) :: y,ymin,ymax
+ real(kind=8) :: wykl=4.0d0
+ if (y.lt.ymin) then
+ gnmr1=(ymin-y)**wykl/wykl
+ else if (y.gt.ymax) then
+ gnmr1=(y-ymax)**wykl/wykl
+ else
+ gnmr1=0.0d0
+ endif
+ return
+ end function gnmr1
+!-----------------------------------------------------------------------------
+ real(kind=8) function gnmr1prim(y,ymin,ymax)
+! implicit none
+ real(kind=8) :: y,ymin,ymax
+ real(kind=8) :: wykl=4.0d0
+ if (y.lt.ymin) then
+ gnmr1prim=-(ymin-y)**(wykl-1)
+ else if (y.gt.ymax) then
+ gnmr1prim=(y-ymax)**(wykl-1)
+ else
+ gnmr1prim=0.0d0
+ endif
+ return
+ end function gnmr1prim
+!----------------------------------------------------------------------------
+ real(kind=8) function rlornmr1(y,ymin,ymax,sigma)
+ real(kind=8) y,ymin,ymax,sigma
+ real(kind=8) wykl /4.0d0/
+ if (y.lt.ymin) then
+ rlornmr1=(ymin-y)**wykl/((ymin-y)**wykl+sigma**wykl)
+ else if (y.gt.ymax) then
+ rlornmr1=(y-ymax)**wykl/((y-ymax)**wykl+sigma**wykl)
+ else
+ rlornmr1=0.0d0
+ endif
+ return
+ end function rlornmr1
+!------------------------------------------------------------------------------
+ real(kind=8) function rlornmr1prim(y,ymin,ymax,sigma)
+ real(kind=8) y,ymin,ymax,sigma
+ real(kind=8) wykl /4.0d0/
+ if (y.lt.ymin) then
+ rlornmr1prim=-(ymin-y)**(wykl-1)*sigma**wykl*wykl/ &
+ ((ymin-y)**wykl+sigma**wykl)**2
+ else if (y.gt.ymax) then
+ rlornmr1prim=(y-ymax)**(wykl-1)*sigma**wykl*wykl/ &
+ ((y-ymax)**wykl+sigma**wykl)**2
+ else
+ rlornmr1prim=0.0d0
+ endif
+ return
+ end function rlornmr1prim
+
+ real(kind=8) function harmonic(y,ymax)
+! implicit none
+ real(kind=8) :: y,ymax
+ real(kind=8) :: wykl=2.0d0
+ harmonic=(y-ymax)**wykl
+ return
+ end function harmonic
+!-----------------------------------------------------------------------------
+ real(kind=8) function harmonicprim(y,ymax)
+ real(kind=8) :: y,ymin,ymax
+ real(kind=8) :: wykl=2.0d0
+ harmonicprim=(y-ymax)*wykl
+ return
+ end function harmonicprim
+!-----------------------------------------------------------------------------
+! gradient_p.F
+!-----------------------------------------------------------------------------
+#ifndef LBFGS
+ subroutine gradient(n,x,nf,g,uiparm,urparm,ufparm)
+
+ use io_base, only:intout,briefout
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.MD'
+! include 'COMMON.IOUNITS'
+ real(kind=8),external :: ufparm
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+ real(kind=8),dimension(6*nres) :: x,g !(maxvar) (maxvar=6*maxres)
+ real(kind=8) :: f,gthetai,gphii,galphai,gomegai
+ integer :: n,nf,ind,ind1,i,k,j
+!
+! This subroutine calculates total internal coordinate gradient.
+! Depending on the number of function evaluations, either whole energy
+! is evaluated beforehand, Cartesian coordinates and their derivatives in
+! internal coordinates are reevaluated or only the cartesian-in-internal
+! coordinate derivatives are evaluated. The subroutine was designed to work
+! with SUMSL.
+!
+!
+ icg=mod(nf,2)+1
+
+!d print *,'grad',nf,icg
+ if (nf-nfl+1) 20,30,40
+ 20 call func(n,x,nf,f,uiparm,urparm,ufparm)
+! write (iout,*) 'grad 20'
+ if (nf.eq.0) return
+ goto 40
+ 30 call var_to_geom(n,x)
+ call chainbuild
+! write (iout,*) 'grad 30'
+!
+! Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+!
+ 40 call cartder
+! write (iout,*) 'grad 40'
+! print *,'GRADIENT: nnt=',nnt,' nct=',nct,' expon=',expon
+!
+! Convert the Cartesian gradient into internal-coordinate gradient.
+!
+ ind=0
+ ind1=0
+ do i=1,nres-2
+ gthetai=0.0D0
+ gphii=0.0D0
+ do j=i+1,nres-1
+ ind=ind+1
+! ind=indmat(i,j)
+! print *,'GRAD: i=',i,' jc=',j,' ind=',ind
+ do k=1,3
+ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
+ enddo
+ do k=1,3
+ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
+ enddo
+ enddo
+ do j=i+1,nres-1
+ ind1=ind1+1
+! ind1=indmat(i,j)
+! print *,'GRAD: i=',i,' jx=',j,' ind1=',ind1
+ do k=1,3
+ gthetai=gthetai+dxdv(k,ind1)*gradx(k,j,icg)
+ gphii=gphii+dxdv(k+3,ind1)*gradx(k,j,icg)
+ enddo
+ enddo
+ if (i.gt.1) g(i-1)=gphii
+ if (n.gt.nphi) g(nphi+i)=gthetai
+ enddo
+ if (n.le.nphi+ntheta) goto 10
+ do i=2,nres-1
+ if (itype(i,1).ne.10) then
+ galphai=0.0D0
+ gomegai=0.0D0
+ do k=1,3
+ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
+ enddo
+ do k=1,3
+ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
+ enddo
+ g(ialph(i,1))=galphai
+ g(ialph(i,1)+nside)=gomegai
+ endif
+ enddo
+!
+! Add the components corresponding to local energy terms.
+!
+ 10 continue
+ do i=1,nvar
+!d write (iout,*) 'i=',i,'g=',g(i),' gloc=',gloc(i,icg)
+ g(i)=g(i)+gloc(i,icg)
+ enddo
+! Uncomment following three lines for diagnostics.
+!d call intout
+!elwrite(iout,*) "in gradient after calling intout"
+!d call briefout(0,0.0d0)
+!d write (iout,'(i3,1pe15.5)') (k,g(k),k=1,n)
+ return
+ end subroutine gradient
+#endif
+!-----------------------------------------------------------------------------
+ subroutine func(n,x,nf,f,uiparm,urparm,ufparm) !from minimize_p.F
+
+ use comm_chu
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.DERIV'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+ integer :: n,nf
+!el integer :: jjj
+!el common /chuju/ jjj
+ real(kind=8) :: energia(0:n_ene)
+ integer :: uiparm(1)
+ real(kind=8) :: urparm(1)
+ real(kind=8) :: f
+ real(kind=8),external :: ufparm
+ real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
+! if (jjj.gt.0) then
+! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+! endif
+ nfl=nf
+ icg=mod(nf,2)+1
+!d print *,'func',nf,nfl,icg
+ call var_to_geom(n,x)
+ call zerograd
+ call chainbuild
+!d write (iout,*) 'ETOTAL called from FUNC'
+ call etotal(energia)
+ call sum_gradient
+ f=energia(0)
+! if (jjj.gt.0) then
+! write (iout,'(10f8.3)') (rad2deg*x(i),i=1,n)
+! write (iout,*) 'f=',etot
+! jjj=0
+! endif
+ return
+ end subroutine func
+!-----------------------------------------------------------------------------
+ subroutine cartgrad
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+ use energy_data
+ use MD_data, only: totT,usampl,eq_time
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.MD'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+!
+ integer :: i,j
+ real(kind=8) :: time00,time01
+
+! This subrouting calculates total Cartesian coordinate gradient.
+! The subroutine chainbuild_cart and energy MUST be called beforehand.
+!
+!#define DEBUG
+#ifdef TIMINGtime01
+ time00=MPI_Wtime()
+#endif
+ icg=1
+ call sum_gradient
+#ifdef TIMING
+#endif
+!#define DEBUG
+!el write (iout,*) "After sum_gradient"
+#ifdef DEBUG
+ write (iout,*) "After sum_gradient"
+ do i=1,nres-1
+ write (iout,*) i," gradc ",(gradc(j,i,icg),j=1,3)
+ write (iout,*) i," gradx ",(gradx(j,i,icg),j=1,3)
+ enddo
+#endif
+!#undef DEBUG
+! If performing constraint dynamics, add the gradients of the constraint energy
+ if(usampl.and.totT.gt.eq_time) then
+ do i=1,nct
+ do j=1,3
+ gradc(j,i,icg)=gradc(j,i,icg)+dudconst(j,i)+duscdiff(j,i)
+ gradx(j,i,icg)=gradx(j,i,icg)+dudxconst(j,i)+duscdiffx(j,i)
+ enddo
+ enddo
+ do i=1,nres-3
+ gloc(i,icg)=gloc(i,icg)+dugamma(i)
+ enddo
+ do i=1,nres-2
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)+dutheta(i)
+ enddo
+ endif
+!elwrite (iout,*) "After sum_gradient"
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
+ call intcartderiv
+!elwrite (iout,*) "After sum_gradient"
+#ifdef TIMING
+ time_intcartderiv=time_intcartderiv+MPI_Wtime()-time01
+#endif
+! call checkintcartgrad
+! write(iout,*) 'calling int_to_cart'
+!#define DEBUG
+#ifdef DEBUG
+ write (iout,*) "gcart, gxcart, gloc before int_to_cart"
+#endif
+ do i=0,nct
+ do j=1,3
+ gcart(j,i)=gradc(j,i,icg)
+ gxcart(j,i)=gradx(j,i,icg)
+! if (i.le.2) print *,"gcart_one",gcart(j,i),gradc(j,i,icg)
+ enddo
+#ifdef DEBUG
+ write (iout,'(i5,2(3f10.5,5x),4f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3),gloc(i,icg),(gloc_sc(j,i,icg),j=1,3)
+#endif
+ enddo
+#ifdef TIMING
+ time01=MPI_Wtime()
+#endif
+! print *,"gcart_two",gxcart(2,2),gradx(2,2,icg)
+ call int_to_cart
+! print *,"gcart_two",gxcart(2,2),gradx(2,2,icg)
+
+#ifdef TIMING
+ time_inttocart=time_inttocart+MPI_Wtime()-time01
+#endif
+#ifdef DEBUG
+ write (iout,*) "gcart and gxcart after int_to_cart"
+ do i=0,nres-1
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),&
+ (gxcart(j,i),j=1,3)
+ enddo
+#endif
+!#undef DEBUG
+#ifdef CARGRAD
+#ifdef DEBUG
+ write (iout,*) "CARGRAD"
+#endif
+! do i=nres,0,-1
+! do j=1,3
+! gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+ ! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+! enddo
+ ! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+ ! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+! enddo
+ ! Correction: dummy residues
+! if (nnt.gt.1) then
+! do j=1,3
+! ! gcart_new(j,nnt)=gcart_new(j,nnt)+gcart_new(j,1)
+! gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+! enddo
+! endif
+! if (nct.lt.nres) then
+! do j=1,3
+! ! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+! gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+! enddo
+! endif
+! call grad_transform
+#endif
+#ifdef TIMING
+ time_cartgrad=time_cartgrad+MPI_Wtime()-time00
+#endif
+!#undef DEBUG
+ return
+ end subroutine cartgrad
+
+#ifdef FIVEDIAG
+ subroutine grad_transform
+ implicit none
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ integer i,j,kk,mnum
+#ifdef DEBUG
+ write (iout,*)"Converting virtual-bond gradient to CA/SC gradient"
+ write (iout,*) "dC/dX gradient"
+ do i=0,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+ & (gxcart(j,i),j=1,3)
+ enddo
+#endif
+ do i=nres,1,-1
+ do j=1,3
+ gcart(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+! gcart_new(j,i)=-gcart(j,i)+gcart(j,i-1)-gxcart(j,i)
+ enddo
+! write (iout,'(i5,3f10.5,5x,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3), &
+! (gcart_new(j,i),j=1,3),(gxcart(j,i),j=1,3)
+ enddo
+! Correction: dummy residues
+ do i=2,nres
+ mnum=molnum(i)
+ if (itype(i-1,mnum).eq.ntyp1_molec(mnum) .and.&
+ itype(i,mnum).ne.ntyp1_molec(mnum)) then
+ gcart(:,i)=gcart(:,i)+gcart(:,i-1)
+ else if (itype(i-1,mnum).ne.ntyp1_molec(mnum).and.&
+ itype(i,mnum).eq.ntyp1_molec(mnum)) then
+ gcart(:,i-1)=gcart(:,i-1)+gcart(:,i)
+ endif
+ enddo
+! if (nnt.gt.1) then
+! do j=1,3
+! gcart(j,nnt)=gcart(j,nnt)+gcart(j,1)
+! enddo
+! endif
+! if (nct.lt.nres) then
+! do j=1,3
+!! gcart_new(j,nct)=gcart_new(j,nct)+gcart_new(j,nres)
+! gcart(j,nct)=gcart(j,nct)+gcart(j,nres)
+! enddo
+! endif
+#ifdef DEBUG
+ write (iout,*) "CA/SC gradient"
+ do i=1,nres
+ write (iout,'(i5,3f10.5,5x,3f10.5)') i,(gcart(j,i),j=1,3),
+ & (gxcart(j,i),j=1,3)
+ enddo
+#endif
+ return
+ end subroutine grad_transform
+#endif
+
+ !-----------------------------------------------------------------------------
+ subroutine zerograd
+ ! implicit real(kind=8) (a-h,o-z)
+ ! include 'DIMENSIONS'
+ ! include 'COMMON.DERIV'
+ ! include 'COMMON.CHAIN'
+ ! include 'COMMON.VAR'
+ ! include 'COMMON.MD'
+ ! include 'COMMON.SCCOR'
+ !
+ !el local variables
+ integer :: i,j,intertyp,k
+ ! Initialize Cartesian-coordinate gradient
+ !
+ ! if (.not.allocated(gradx)) allocate(gradx(3,nres,2)) !(3,maxres,2)
+ ! if (.not.allocated(gradc)) allocate(gradc(3,nres,2)) !(3,maxres,2)
+
+ ! allocate(gvdwx(3,nres),gvdwc(3,nres),gelc(3,nres),gelc_long(3,nres))
+ ! allocate(gvdwpp(3,nres),gvdwc_scpp(3,nres),gradx_scp(3,nres))
+ ! allocate(gvdwc_scp(3,nres),ghpbx(3,nres),ghpbc(3,nres))
+ ! allocate(gradcorr_long(3,nres))
+ ! allocate(gradcorr5_long(3,nres),gradcorr6_long(3,nres))
+ ! allocate(gcorr6_turn_long(3,nres))
+ ! allocate(gradcorr5(3,nres),gradcorr6(3,nres)) !(3,maxres)
+
+ ! if (.not.allocated(gradcorr)) allocate(gradcorr(3,nres)) !(3,maxres)
+
+ ! allocate(gel_loc(3,nres),gel_loc_long(3,nres),gcorr3_turn(3,nres))
+ ! allocate(gcorr4_turn(3,nres),gcorr6_turn(3,nres))
+
+ ! if (.not.allocated(gradb)) allocate(gradb(3,nres)) !(3,maxres)
+ ! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
+
+ ! allocate(gsccorc(3,nres),gsccorx(3,nres)) !(3,maxres)
+ ! allocate(gscloc(3,nres)) !(3,maxres)
+ ! if (.not.allocated(gsclocx)) allocate(gsclocx(3,nres)) !(3,maxres)
+
+
+
+ ! common /deriv_scloc/
+ ! allocate(dXX_C1tab(3,nres),dYY_C1tab(3,nres),dZZ_C1tab(3,nres))
+ ! allocate(dXX_Ctab(3,nres),dYY_Ctab(3,nres),dZZ_Ctab(3,nres))
+ ! allocate(dXX_XYZtab(3,nres),dYY_XYZtab(3,nres),dZZ_XYZtab(3,nres)) !(3,maxres)
+ ! common /mpgrad/
+ ! allocate(jgrad_start(nres),jgrad_end(nres)) !(maxres)
+
+
+
+ ! gradc(j,i,icg)=0.0d0
+ ! gradx(j,i,icg)=0.0d0
+
+ ! allocate(gloc_sc(3,nres,10)) !(3,0:maxres2,10)maxres2=2*maxres
+ !elwrite(iout,*) "icg",icg
+ do i=-1,nres
+ do j=1,3
+ gvdwx(j,i)=0.0D0
+ gradx_scp(j,i)=0.0D0
+ gvdwc(j,i)=0.0D0
+ gvdwc_scp(j,i)=0.0D0
+ gvdwc_scpp(j,i)=0.0d0
+ gelc(j,i)=0.0D0
+ gelc_long(j,i)=0.0D0
+ gradb(j,i)=0.0d0
+ gradbx(j,i)=0.0d0
+ gvdwpp(j,i)=0.0d0
+ gel_loc(j,i)=0.0d0
+ gel_loc_long(j,i)=0.0d0
+ ghpbc(j,i)=0.0D0
+ ghpbx(j,i)=0.0D0
+ gcorr3_turn(j,i)=0.0d0
+ gcorr4_turn(j,i)=0.0d0
+ gradcorr(j,i)=0.0d0
+ gradcorr_long(j,i)=0.0d0
+ gradcorr5_long(j,i)=0.0d0
+ gradcorr6_long(j,i)=0.0d0
+ gcorr6_turn_long(j,i)=0.0d0
+ gradcorr5(j,i)=0.0d0
+ gradcorr6(j,i)=0.0d0
+ gcorr6_turn(j,i)=0.0d0
+ gsccorc(j,i)=0.0d0
+ gsccorx(j,i)=0.0d0
+ gradc(j,i,icg)=0.0d0
+ gradx(j,i,icg)=0.0d0
+ gscloc(j,i)=0.0d0
+ gsclocx(j,i)=0.0d0
+ gliptran(j,i)=0.0d0
+ gliptranx(j,i)=0.0d0
+ gliptranc(j,i)=0.0d0
+ gshieldx(j,i)=0.0d0
+ gshieldc(j,i)=0.0d0
+ gshieldc_loc(j,i)=0.0d0
+ gshieldx_ec(j,i)=0.0d0
+ gshieldc_ec(j,i)=0.0d0
+ gshieldc_loc_ec(j,i)=0.0d0
+ gshieldx_t3(j,i)=0.0d0
+ gshieldc_t3(j,i)=0.0d0
+ gshieldc_loc_t3(j,i)=0.0d0
+ gshieldx_t4(j,i)=0.0d0
+ gshieldc_t4(j,i)=0.0d0
+ gshieldc_loc_t4(j,i)=0.0d0
+ gshieldx_ll(j,i)=0.0d0
+ gshieldc_ll(j,i)=0.0d0
+ gshieldc_loc_ll(j,i)=0.0d0
+ gg_tube(j,i)=0.0d0
+ gg_tube_sc(j,i)=0.0d0
+ gradafm(j,i)=0.0d0
+ gradb_nucl(j,i)=0.0d0
+ gradbx_nucl(j,i)=0.0d0
+ gvdwpp_nucl(j,i)=0.0d0
+ gvdwpp(j,i)=0.0d0
+ gelpp(j,i)=0.0d0
+ gvdwpsb(j,i)=0.0d0
+ gvdwpsb1(j,i)=0.0d0
+ gvdwsbc(j,i)=0.0d0
+ gvdwsbx(j,i)=0.0d0
+ gelsbc(j,i)=0.0d0
+ gradcorr_nucl(j,i)=0.0d0
+ gradcorr3_nucl(j,i)=0.0d0
+ gradxorr_nucl(j,i)=0.0d0
+ gradxorr3_nucl(j,i)=0.0d0
+ gelsbx(j,i)=0.0d0
+ gsbloc(j,i)=0.0d0
+ gsblocx(j,i)=0.0d0
+ gradpepcat(j,i)=0.0d0
+ gradpepcatx(j,i)=0.0d0
+ gradcatcat(j,i)=0.0d0
+ gvdwx_scbase(j,i)=0.0d0
+ gvdwc_scbase(j,i)=0.0d0
+ gvdwx_pepbase(j,i)=0.0d0
+ gvdwc_pepbase(j,i)=0.0d0
+ gvdwx_scpho(j,i)=0.0d0
+ gvdwc_scpho(j,i)=0.0d0
+ gvdwc_peppho(j,i)=0.0d0
+ gradnuclcatx(j,i)=0.0d0
+ gradnuclcat(j,i)=0.0d0
+ gradlipbond(j,i)=0.0d0
+ gradlipang(j,i)=0.0d0
+ gradliplj(j,i)=0.0d0
+ gradlipelec(j,i)=0.0d0
+ gradcattranc(j,i)=0.0d0
+ gradcattranx(j,i)=0.0d0
+ gradcatangx(j,i)=0.0d0
+ gradcatangc(j,i)=0.0d0
+ gradpepmart(j,i)=0.0d0
+ gradpepmartx(j,i)=0.0d0
+ duscdiff(j,i)=0.0d0
+ duscdiffx(j,i)=0.0d0
+ enddo
+ enddo
+ do i=0,nres
+ do j=1,3
+ do intertyp=1,3
+ gloc_sc(intertyp,i,icg)=0.0d0
+ enddo
+ enddo
+ enddo
+ do i=1,nres
+ do j=1,maxcontsshi
+ shield_list(j,i)=0
+ do k=1,3
+ !C print *,i,j,k
+ grad_shield_side(k,j,i)=0.0d0
+ grad_shield_loc(k,j,i)=0.0d0
+ enddo
+ enddo
+ ishield_list(i)=0
+ enddo
+
+ !
+ ! Initialize the gradient of local energy terms.
+ !
+ ! allocate(gloc(4*nres,2)) !!(maxvar,2)(maxvar=6*maxres)
+ ! if (.not.allocated(gel_loc_loc)) allocate(gel_loc_loc(nres)) !(maxvar)(maxvar=6*maxres)
+ ! if (.not.allocated(gcorr_loc)) allocate(gcorr_loc(nres)) !(maxvar)(maxvar=6*maxres)
+ ! allocate(g_corr5_loc(nres),g_corr6_loc(nres)) !(maxvar)(maxvar=6*maxres)
+ ! allocate(gel_loc_turn3(nres))
+ ! allocate(gel_loc_turn4(nres),gel_loc_turn6(nres)) !(maxvar)(maxvar=6*maxres)
+ ! allocate(gsccor_loc(nres)) !(maxres)
+
+ do i=1,4*nres
+ gloc(i,icg)=0.0D0
+ enddo
+ do i=1,nres
+ gel_loc_loc(i)=0.0d0
+ gcorr_loc(i)=0.0d0
+ g_corr5_loc(i)=0.0d0
+ g_corr6_loc(i)=0.0d0
+ gel_loc_turn3(i)=0.0d0
+ gel_loc_turn4(i)=0.0d0
+ gel_loc_turn6(i)=0.0d0
+ gsccor_loc(i)=0.0d0
+ enddo
+ ! initialize gcart and gxcart
+ ! allocate(gcart(3,0:nres),gxcart(3,0:nres)) !(3,0:MAXRES)
+ do i=0,nres
+ do j=1,3
+ gcart(j,i)=0.0d0
+ gxcart(j,i)=0.0d0
+ enddo
+ enddo
+ return
+ end subroutine zerograd
+ !-----------------------------------------------------------------------------
+ real(kind=8) function fdum()
+ fdum=0.0D0
+ return
+ end function fdum
+ !-----------------------------------------------------------------------------
+ ! intcartderiv.F
+ !-----------------------------------------------------------------------------
+ subroutine intcartderiv
+ ! implicit real(kind=8) (a-h,o-z)
+ ! include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ ! include 'COMMON.SETUP'
+ ! include 'COMMON.CHAIN'
+ ! include 'COMMON.VAR'
+ ! include 'COMMON.GEO'
+ ! include 'COMMON.INTERACT'
+ ! include 'COMMON.DERIV'
+ ! include 'COMMON.IOUNITS'
+ ! include 'COMMON.LOCAL'
+ ! include 'COMMON.SCCOR'
+ real(kind=8) :: pi4,pi34
+ real(kind=8),dimension(3,2,nres) :: dcostheta ! (3,2,maxres)
+ real(kind=8),dimension(3,3,nres) :: dcosphi,dsinphi,dcosalpha,&
+ dcosomega,dsinomega !(3,3,maxres)
+ real(kind=8),dimension(3) :: vo1,vo2,vo3,dummy,vp1,vp2,vp3,vpp1,n
+
+ integer :: i,j,k
+ real(kind=8) :: cost,sint,cost1,sint1,cost2,sint2,sing,cosg,scalp,&
+ fac0,fac1,fac2,fac3,fac4,fac5,fac6,ctgt,ctgt1,cosg_inv,&
+ fac7,fac8,fac9,scala1,scala2,cosa,sina,sino,fac15,fac16,&
+ fac17,coso_inv,fac10,fac11,fac12,fac13,fac14,IERROR
+ integer :: nres2
+ nres2=2*nres
+
+ !el from module energy-------------
+ !el allocate(dcostau(3,3,3,itau_start:itau_end)) !(3,3,3,maxres2)maxres2=2*maxres
+ !el allocate(dsintau(3,3,3,itau_start:itau_end))
+ !el allocate(dtauangle(3,3,3,itau_start:itau_end))
+
+ !el allocate(dcostau(3,3,3,0:nres2)) !(3,3,3,maxres2)maxres2=2*maxres
+ !el allocate(dsintau(3,3,3,0:nres2))
+ !el allocate(dtauangle(3,3,3,0:nres2))
+ !el allocate(domicron(3,2,2,0:nres2))
+ !el allocate(dcosomicron(3,2,2,0:nres2))
+
+
+
+#if defined(MPI) && defined(PARINTDER)
+ if (nfgtasks.gt.1 .and. me.eq.king) &
+ call MPI_Bcast(8,1,MPI_INTEGER,king,FG_COMM,IERROR)
+#endif
+ pi4 = 0.5d0*pipol
+ pi34 = 3*pi4
+
+ ! allocate(dtheta(3,2,nres)) !(3,2,maxres)
+ ! allocate(dphi(3,3,nres),dalpha(3,3,nres),domega(3,3,nres)) !(3,3,maxres)
+
+ ! write (iout,*) "iphi1_start",iphi1_start," iphi1_end",iphi1_end
+ do i=1,nres
+ do j=1,3
+ dtheta(j,1,i)=0.0d0
+ dtheta(j,2,i)=0.0d0
+ dphi(j,1,i)=0.0d0
+ dphi(j,2,i)=0.0d0
+ dphi(j,3,i)=0.0d0
+ dcosomicron(j,1,1,i)=0.0d0
+ dcosomicron(j,1,2,i)=0.0d0
+ dcosomicron(j,2,1,i)=0.0d0
+ dcosomicron(j,2,2,i)=0.0d0
+ enddo
+ enddo
+ ! Derivatives of theta's
+#if defined(MPI) && defined(PARINTDER)
+ ! We need dtheta(:,:,i-1) to compute dphi(:,:,i)
+ do i=max0(ithet_start-1,3),ithet_end
+#else
+ do i=3,nres
+#endif
+ cost=dcos(theta(i))
+ sint=sqrt(1-cost*cost)
+ do j=1,3
+ dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/&
+ vbld(i-1)
+ if (((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))) &
+ dtheta(j,1,i)=-dcostheta(j,1,i)/sint
+ dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/&
+ vbld(i)
+ if ((itype(i-1,1).ne.ntyp1).and.(sint.ne.0.0d0))&
+ dtheta(j,2,i)=-dcostheta(j,2,i)/sint
+ enddo
+ enddo
+#if defined(MPI) && defined(PARINTDER)
+ ! We need dtheta(:,:,i-1) to compute dphi(:,:,i)
+ do i=max0(ithet_start-1,3),ithet_end
+#else
+ do i=3,nres
+#endif
+ if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1).and.molnum(i).lt.4) then
+ cost1=dcos(omicron(1,i))
+ sint1=sqrt(1-cost1*cost1)
+ cost2=dcos(omicron(2,i))
+ sint2=sqrt(1-cost2*cost2)
+ do j=1,3
+ !C Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1)
+ dcosomicron(j,1,1,i)=-(dc_norm(j,i-1+nres)+ &
+ cost1*dc_norm(j,i-2))/ &
+ vbld(i-1)
+ domicron(j,1,1,i)=-1.0/sint1*dcosomicron(j,1,1,i)
+ dcosomicron(j,1,2,i)=-(dc_norm(j,i-2) &
+ +cost1*(dc_norm(j,i-1+nres)))/ &
+ vbld(i-1+nres)
+ domicron(j,1,2,i)=-1.0/sint1*dcosomicron(j,1,2,i)
+ !C Calculate derivative over second omicron Sci-1,Cai-1 Cai
+ !C Looks messy but better than if in loop
+ dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres) &
+ +cost2*dc_norm(j,i-1))/ &
+ vbld(i)
+ domicron(j,2,1,i)=-1.0/sint2*dcosomicron(j,2,1,i)
+ dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) &
+ +cost2*(-dc_norm(j,i-1+nres)))/ &
+ vbld(i-1+nres)
+ ! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres)
+ domicron(j,2,2,i)=-1.0/sint2*dcosomicron(j,2,2,i)
+ enddo
+ endif
+ enddo
+ !elwrite(iout,*) "after vbld write"
+ ! Derivatives of phi:
+ ! If phi is 0 or 180 degrees, then the formulas
+ ! have to be derived by power series expansion of the
+ ! conventional formulas around 0 and 180.
+#ifdef PARINTDER
+ do i=iphi1_start,iphi1_end
+#else
+ do i=4,nres
+#endif
+ ! if (itype(i-1,1).eq.21 .or. itype(i-2,1).eq.21 ) cycle
+ ! the conventional case
+ sint=dsin(theta(i))
+ sint1=dsin(theta(i-1))
+ sing=dsin(phi(i))
+ cost=dcos(theta(i))
+ cost1=dcos(theta(i-1))
+ cosg=dcos(phi(i))
+ scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1))
+ if ((sint*sint1).eq.0.0d0) then
+ fac0=0.0d0
+ else
+ fac0=1.0d0/(sint1*sint)
+ endif
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ if (sint1.ne.0.0d0) then
+ fac3=cosg*cost1/(sint1*sint1)
+ else
+ fac3=0.0d0
+ endif
+ if (sint.ne.0.0d0) then
+ fac4=cosg*cost/(sint*sint)
+ else
+ fac4=0.0d0
+ endif
+ ! Obtaining the gamma derivatives from sine derivative
+ if (phi(i).gt.-pi4.and.phi(i).le.pi4.or. &
+ phi(i).gt.pi34.and.phi(i).le.pi.or. &
+ phi(i).ge.-pi.and.phi(i).le.-pi34) then
+ call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-1),vp2)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
+ do j=1,3
+ if (sint.ne.0.0d0) then
+ ctgt=cost/sint
+ else
+ ctgt=0.0d0
+ endif
+ if (sint1.ne.0.0d0) then
+ ctgt1=cost1/sint1
+ else
+ ctgt1=0.0d0
+ endif
+ cosg_inv=1.0d0/cosg
+! if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
+ dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
+ -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
+ dphi(j,1,i)=cosg_inv*dsinphi(j,1,i)
+ dsinphi(j,2,i)= &
+ -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*dtheta(j,1,i)) &
+ -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ dphi(j,2,i)=cosg_inv*dsinphi(j,2,i)
+ dsinphi(j,3,i)=-sing*ctgt*dtheta(j,2,i) &
+ +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
+ ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dphi(j,3,i)=cosg_inv*dsinphi(j,3,i)
+! endif
+! write(iout,*) "just after,close to pi",dphi(j,3,i),&
+! sing*(ctgt1*dtheta(j,2,i-1)),ctgt*dtheta(j,1,i), &
+! (fac0*vp2(j)+sing*dc_norm(j,i-2)),vbld_inv(i-1)
+
+ ! Bug fixed 3/24/05 (AL)
+ enddo
+ ! Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+! if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then
+ dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
+ dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
+ dc_norm(j,i-3))/vbld(i-2)
+ dphi(j,1,i)=-1.0/sing*dcosphi(j,1,i)
+ dcosphi(j,2,i)=fac1*dcostheta(j,2,i-1)+fac2* &
+ dcostheta(j,1,i)+fac3*dcostheta(j,2,i-1)+fac4* &
+ dcostheta(j,1,i)
+ dphi(j,2,i)=-1.0/sing*dcosphi(j,2,i)
+ dcosphi(j,3,i)=fac2*dcostheta(j,2,i)+fac4* &
+ dcostheta(j,2,i)-fac0*(dc_norm(j,i-3)-scalp* &
+ dc_norm(j,i-1))/vbld(i)
+ dphi(j,3,i)=-1.0/sing*dcosphi(j,3,i)
+!#define DEBUG
+#ifdef DEBUG
+ write(iout,*) "just after",dphi(j,3,i),sing,dcosphi(j,3,i)
+#endif
+!#undef DEBUG
+! endif
+ enddo
+ endif
+ enddo
+ !alculate derivative of Tauangle
+#ifdef PARINTDER
+ do i=itau_start,itau_end
+#else
+ do i=3,nres
+ !elwrite(iout,*) " vecpr",i,nres
+#endif
+ if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
+ ! if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10).or.
+ ! & (itype(i-1,1).eq.ntyp1).or.(itype(i,1).eq.ntyp1)) cycle
+ !c dtauangle(j,intertyp,dervityp,residue number)
+ !c INTERTYP=1 SC...Ca...Ca..Ca
+ ! the conventional case
+ sint=dsin(theta(i))
+ sint1=dsin(omicron(2,i-1))
+ sing=dsin(tauangle(1,i))
+ cost=dcos(theta(i))
+ cost1=dcos(omicron(2,i-1))
+ cosg=dcos(tauangle(1,i))
+ !elwrite(iout,*) " vecpr5",i,nres
+ do j=1,3
+ !elwrite(iout,*) " vecpreee",i,nres,j,i-2+nres
+ !elwrite(iout,*) " vecpr5",dc_norm2(1,1)
+ dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+ ! write(iout,*) dc_norm2(j,i-2+nres),"dcnorm"
+ enddo
+ scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
+ ! write(iout,*) "faki",fac0,fac1,fac2,fac3,fac
+ if ((sint*sint1).eq.0.0d0) then
+ fac0=0.0d0
+ else
+ fac0=1.0d0/(sint1*sint)
+ endif
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ if (sint1.ne.0.0d0) then
+ fac3=cosg*cost1/(sint1*sint1)
+ else
+ fac3=0.0d0
+ endif
+ if (sint.ne.0.0d0) then
+ fac4=cosg*cost/(sint*sint)
+ else
+ fac4=0.0d0
+ endif
+
+ ! Obtaining the gamma derivatives from sine derivative
+ if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or. &
+ tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or. &
+ tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
+ call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,2,i-1) &
+ -(fac0*vp1(j)+sing*(dc_norm2(j,i-2+nres))) &
+ *vbld_inv(i-2+nres)
+ dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
+ dsintau(j,1,2,i)= &
+ -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*dtheta(j,1,i)) &
+ -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ ! write(iout,*) "dsintau", dsintau(j,1,2,i)
+ dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
+ ! Bug fixed 3/24/05 (AL)
+ dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i) &
+ +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
+ ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dtauangle(j,1,3,i)=cosg_inv*dsintau(j,1,3,i)
+ enddo
+ ! Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcostau(j,1,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3* &
+ dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1)-scalp* &
+ (dc_norm2(j,i-2+nres)))/vbld(i-2+nres)
+ dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
+ dcostau(j,1,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2* &
+ dcostheta(j,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4* &
+ dcostheta(j,1,i)
+ dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
+ dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4* &
+ dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp* &
+ dc_norm(j,i-1))/vbld(i)
+ dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
+ ! write (iout,*) "else",i
+ enddo
+ endif
+ ! do k=1,3
+ ! write(iout,*) "tu",i,k,(dtauangle(j,1,k,i),j=1,3)
+ ! enddo
+ enddo
+ !C Second case Ca...Ca...Ca...SC
+#ifdef PARINTDER
+ do i=itau_start,itau_end
+#else
+ do i=4,nres
+#endif
+ if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+ (itype(i-2,1).eq.ntyp1).or.(itype(i-3,1).eq.ntyp1)) cycle
+ ! the conventional case
+ sint=dsin(omicron(1,i))
+ sint1=dsin(theta(i-1))
+ sing=dsin(tauangle(2,i))
+ cost=dcos(omicron(1,i))
+ cost1=dcos(theta(i-1))
+ cosg=dcos(tauangle(2,i))
+ ! do j=1,3
+ ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+ ! enddo
+ scalp=scalar(dc_norm(1,i-3),dc_norm(1,i-1+nres))
+ if ((sint*sint1).eq.0.0d0) then
+ fac0=0.0d0
+ else
+ fac0=1.0d0/(sint1*sint)
+ endif
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ if (sint1.ne.0.0d0) then
+ fac3=cosg*cost1/(sint1*sint1)
+ else
+ fac3=0.0d0
+ endif
+ if (sint.ne.0.0d0) then
+ fac4=cosg*cost/(sint*sint)
+ else
+ fac4=0.0d0
+ endif
+ ! Obtaining the gamma derivatives from sine derivative
+ if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or. &
+ tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or. &
+ tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
+ call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-1+nres),vp2)
+ call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1) &
+ +(fac0*vp1(j)-sing*dc_norm(j,i-3))*vbld_inv(i-2)
+ ! write(iout,*) i,j,dsintau(j,2,1,i),sing*ctgt1*dtheta(j,1,i-1),
+ ! &fac0*vp1(j),sing*dc_norm(j,i-3),vbld_inv(i-2),"dsintau(2,1)"
+ dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
+ dsintau(j,2,2,i)= &
+ -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i)) &
+ -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ ! write(iout,*) "sprawdzenie",i,j,sing*ctgt1*dtheta(j,2,i-1),
+ ! & sing*ctgt*domicron(j,1,2,i),
+ ! & (fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ dtauangle(j,2,2,i)=cosg_inv*dsintau(j,2,2,i)
+ ! Bug fixed 3/24/05 (AL)
+ dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i) &
+ +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres))*vbld_inv(i-1+nres)
+ ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dtauangle(j,2,3,i)=cosg_inv*dsintau(j,2,3,i)
+ enddo
+ ! Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3* &
+ dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp* &
+ dc_norm(j,i-3))/vbld(i-2)
+ dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
+ dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2* &
+ dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4* &
+ dcosomicron(j,1,1,i)
+ dtauangle(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
+ dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4* &
+ dcosomicron(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp* &
+ dc_norm(j,i-1+nres))/vbld(i-1+nres)
+ dtauangle(j,2,3,i)=-1/sing*dcostau(j,2,3,i)
+ ! write(iout,*) i,j,"else", dtauangle(j,2,3,i)
+ enddo
+ endif
+ enddo
+
+ !CC third case SC...Ca...Ca...SC
+#ifdef PARINTDER
+
+ do i=itau_start,itau_end
+#else
+ do i=3,nres
+#endif
+ ! the conventional case
+ if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. &
+ (itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle
+ sint=dsin(omicron(1,i))
+ sint1=dsin(omicron(2,i-1))
+ sing=dsin(tauangle(3,i))
+ cost=dcos(omicron(1,i))
+ cost1=dcos(omicron(2,i-1))
+ cosg=dcos(tauangle(3,i))
+ do j=1,3
+ dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+ ! dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+ enddo
+ scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres))
+ if ((sint*sint1).eq.0.0d0) then
+ fac0=0.0d0
+ else
+ fac0=1.0d0/(sint1*sint)
+ endif
+ fac1=cost*fac0
+ fac2=cost1*fac0
+ if (sint1.ne.0.0d0) then
+ fac3=cosg*cost1/(sint1*sint1)
+ else
+ fac3=0.0d0
+ endif
+ if (sint.ne.0.0d0) then
+ fac4=cosg*cost/(sint*sint)
+ else
+ fac4=0.0d0
+ endif
+ ! Obtaining the gamma derivatives from sine derivative
+ if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or. &
+ tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or. &
+ tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
+ call vecpr(dc_norm(1,i-1+nres),dc_norm(1,i-2),vp1)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1+nres),vp2)
+ call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+ do j=1,3
+ ctgt=cost/sint
+ ctgt1=cost1/sint1
+ cosg_inv=1.0d0/cosg
+ dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,2,i-1) &
+ -(fac0*vp1(j)-sing*dc_norm(j,i-2+nres)) &
+ *vbld_inv(i-2+nres)
+ dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
+ dsintau(j,3,2,i)= &
+ -sing*(ctgt1*domicron(j,2,1,i-1)+ctgt*domicron(j,1,1,i)) &
+ -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+ dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
+ ! Bug fixed 3/24/05 (AL)
+ dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i) &
+ +(fac0*vp3(j)-sing*dc_norm(j,i-1+nres)) &
+ *vbld_inv(i-1+nres)
+ ! & +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+ dtauangle(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
+ enddo
+ ! Obtaining the gamma derivatives from cosine derivative
+ else
+ do j=1,3
+ dcostau(j,3,1,i)=fac1*dcosomicron(j,2,2,i-1)+fac3* &
+ dcosomicron(j,2,2,i-1)-fac0*(dc_norm(j,i-1+nres)-scalp* &
+ dc_norm2(j,i-2+nres))/vbld(i-2+nres)
+ dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
+ dcostau(j,3,2,i)=fac1*dcosomicron(j,2,1,i-1)+fac2* &
+ dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,1,i-1)+fac4* &
+ dcosomicron(j,1,1,i)
+ dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
+ dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4* &
+ dcosomicron(j,1,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp* &
+ dc_norm(j,i-1+nres))/vbld(i-1+nres)
+ dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
+ ! write(iout,*) "else",i
+ enddo
+ endif
+ enddo
+
+#ifdef CRYST_SC
+ ! Derivatives of side-chain angles alpha and omega
+#if defined(MPI) && defined(PARINTDER)
+ do i=ibond_start,ibond_end
+#else
+ do i=2,nres-1
+#endif
+ if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
+ fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1))))
+ fac6=fac5/vbld(i)
+ fac7=fac5*fac5
+ fac8=fac5/vbld(i+1)
+ fac9=fac5/vbld(i+nres)
+ scala1=scalar(dc_norm(1,i-1),dc_norm(1,i+nres))
+ scala2=scalar(dc_norm(1,i),dc_norm(1,i+nres))
+ cosa=dsqrt(0.5d0/(1.0d0+dcos(theta(i+1))))* &
+ (scalar(dC_norm(1,i),dC_norm(1,i+nres)) &
+ -scalar(dC_norm(1,i-1),dC_norm(1,i+nres)))
+ sina=sqrt(1-cosa*cosa)
+ sino=dsin(omeg(i))
+ ! write (iout,*) "i",i," cosa",cosa," sina",sina," sino",sino
+ do j=1,3
+ dcosalpha(j,1,i)=fac6*(scala1*dc_norm(j,i-1)- &
+ dc_norm(j,i+nres))-cosa*fac7*dcostheta(j,1,i+1)
+ dalpha(j,1,i)=-1/sina*dcosalpha(j,1,i)
+ dcosalpha(j,2,i)=fac8*(dc_norm(j,i+nres)- &
+ scala2*dc_norm(j,i))-cosa*fac7*dcostheta(j,2,i+1)
+ dalpha(j,2,i)=-1/sina*dcosalpha(j,2,i)
+ dcosalpha(j,3,i)=(fac9*(dc_norm(j,i)- &
+ dc_norm(j,i-1))-(cosa*dc_norm(j,i+nres))/ &
+ vbld(i+nres))
+ dalpha(j,3,i)=-1/sina*dcosalpha(j,3,i)
+ enddo
+ ! obtaining the derivatives of omega from sines
+ if(omeg(i).gt.-pi4.and.omeg(i).le.pi4.or. &
+ omeg(i).gt.pi34.and.omeg(i).le.pi.or. &
+ omeg(i).gt.-pi.and.omeg(i).le.-pi34) then
+ fac15=dcos(theta(i+1))/(dsin(theta(i+1))* &
+ dsin(theta(i+1)))
+ fac16=dcos(alph(i))/(dsin(alph(i))*dsin(alph(i)))
+ fac17=1.0d0/(dsin(theta(i+1))*dsin(alph(i)))
+ call vecpr(dc_norm(1,i+nres),dc_norm(1,i),vo1)
+ call vecpr(dc_norm(1,i+nres),dc_norm(1,i-1),vo2)
+ call vecpr(dc_norm(1,i),dc_norm(1,i-1),vo3)
+ coso_inv=1.0d0/dcos(omeg(i))
+ do j=1,3
+ dsinomega(j,1,i)=sino*(fac15*dcostheta(j,1,i+1) &
+ +fac16*dcosalpha(j,1,i))-fac17/vbld(i)*vo1(j)- &
+ (sino*dc_norm(j,i-1))/vbld(i)
+ domega(j,1,i)=coso_inv*dsinomega(j,1,i)
+ dsinomega(j,2,i)=sino*(fac15*dcostheta(j,2,i+1) &
+ +fac16*dcosalpha(j,2,i))+fac17/vbld(i+1)*vo2(j) &
+ -sino*dc_norm(j,i)/vbld(i+1)
+ domega(j,2,i)=coso_inv*dsinomega(j,2,i)
+ dsinomega(j,3,i)=sino*fac16*dcosalpha(j,3,i)- &
+ fac17/vbld(i+nres)*vo3(j)-sino*dc_norm(j,i+nres)/ &
+ vbld(i+nres)
+ domega(j,3,i)=coso_inv*dsinomega(j,3,i)
+ enddo
+ else
+ ! obtaining the derivatives of omega from cosines
+ fac10=sqrt(0.5d0*(1-dcos(theta(i+1))))
+ fac11=sqrt(0.5d0*(1+dcos(theta(i+1))))
+ fac12=fac10*sina
+ fac13=fac12*fac12
+ fac14=sina*sina
+ do j=1,3
+ dcosomega(j,1,i)=(-(0.25d0*cosa/fac11* &
+ dcostheta(j,1,i+1)+fac11*dcosalpha(j,1,i))*fac12+ &
+ (0.25d0/fac10*sina*dcostheta(j,1,i+1)+cosa/sina* &
+ fac10*dcosalpha(j,1,i))*(scala2-fac11*cosa))/fac13
+ domega(j,1,i)=-1/sino*dcosomega(j,1,i)
+ dcosomega(j,2,i)=(((dc_norm(j,i+nres)-scala2* &
+ dc_norm(j,i))/vbld(i+1)-0.25d0*cosa/fac11* &
+ dcostheta(j,2,i+1)-fac11*dcosalpha(j,2,i))*fac12+ &
+ (scala2-fac11*cosa)*(0.25d0*sina/fac10* &
+ dcostheta(j,2,i+1)+fac10*cosa/sina*dcosalpha(j,2,i)))/fac13
+ domega(j,2,i)=-1/sino*dcosomega(j,2,i)
+ dcosomega(j,3,i)=1/fac10*((1/vbld(i+nres)*(dc_norm(j,i)- &
+ scala2*dc_norm(j,i+nres))-fac11*dcosalpha(j,3,i))*sina+ &
+ (scala2-fac11*cosa)*(cosa/sina*dcosalpha(j,3,i)))/fac14
+ domega(j,3,i)=-1/sino*dcosomega(j,3,i)
+ enddo
+ endif
+ else
+ do j=1,3
+ do k=1,3
+ dalpha(k,j,i)=0.0d0
+ domega(k,j,i)=0.0d0
+ enddo
+ enddo
+ endif
+ enddo
+#endif
+#if defined(MPI) && defined(PARINTDER)
+ if (nfgtasks.gt.1) then
+#ifdef DEBUG
+ !d write (iout,*) "Gather dtheta"
+ !d call flush(iout)
+ write (iout,*) "dtheta before gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),k=1,3),j=1,2)
+ enddo
+#endif
+ call MPI_Gatherv(dtheta(1,1,ithet_start),ithet_count(fg_rank),&
+ MPI_THET,dtheta(1,1,1),ithet_count(0),ithet_displ(0),MPI_THET,&
+ king,FG_COMM,IERROR)
+!#define DEBUG
+#ifdef DEBUG
+ !d write (iout,*) "Gather dphi"
+ !d call flush(iout)
+ write (iout,*) "dphi before gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),k=1,3),j=1,3)
+ enddo
+#endif
+!#undef DEBUG
+ call MPI_Gatherv(dphi(1,1,iphi1_start),iphi1_count(fg_rank),&
+ MPI_GAM,dphi(1,1,1),iphi1_count(0),iphi1_displ(0),MPI_GAM,&
+ king,FG_COMM,IERROR)
+ !d write (iout,*) "Gather dalpha"
+ !d call flush(iout)
+#ifdef CRYST_SC
+ call MPI_Gatherv(dalpha(1,1,ibond_start),ibond_count(fg_rank),&
+ MPI_GAM,dalpha(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,&
+ king,FG_COMM,IERROR)
+ !d write (iout,*) "Gather domega"
+ !d call flush(iout)
+ call MPI_Gatherv(domega(1,1,ibond_start),ibond_count(fg_rank),&
+ MPI_GAM,domega(1,1,1),ibond_count(0),ibond_displ(0),MPI_GAM,&
+ king,FG_COMM,IERROR)
+#endif
+ endif
+#endif
+!#define DEBUG
+#ifdef DEBUG
+ write (iout,*) "dtheta after gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dtheta(j,k,i),j=1,3),k=1,2)
+ enddo
+ write (iout,*) "dphi after gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dphi(j,k,i),j=1,3),k=1,3)
+ enddo
+ write (iout,*) "dalpha after gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((dalpha(j,k,i),j=1,3),k=1,3)
+ enddo
+ write (iout,*) "domega after gather"
+ do i=1,nres
+ write (iout,'(i3,3(3f8.5,3x))') i,((domega(j,k,i),j=1,3),k=1,3)
+ enddo
+#endif
+!#undef DEBUG
+ return
+ end subroutine intcartderiv
+ !-----------------------------------------------------------------------------
+ subroutine checkintcartgrad
+ ! implicit real(kind=8) (a-h,o-z)
+ ! include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+ ! include 'COMMON.CHAIN'
+ ! include 'COMMON.VAR'
+ ! include 'COMMON.GEO'
+ ! include 'COMMON.INTERACT'
+ ! include 'COMMON.DERIV'
+ ! include 'COMMON.IOUNITS'
+ ! include 'COMMON.SETUP'
+ real(kind=8),dimension(3,2,nres) :: dthetanum !(3,2,maxres)
+ real(kind=8),dimension(3,3,nres) :: dphinum,dalphanum,domeganum !(3,3,maxres)
+ real(kind=8),dimension(nres) :: theta_s,phi_s,alph_s,omeg_s !(maxres)
+ real(kind=8),dimension(3) :: dc_norm_s
+ real(kind=8) :: aincr=1.0d-5
+ integer :: i,j
+ real(kind=8) :: dcji
+ do i=1,nres
+ phi_s(i)=phi(i)
+ theta_s(i)=theta(i)
+ alph_s(i)=alph(i)
+ omeg_s(i)=omeg(i)
+ enddo
+ ! Check theta gradient
+ write (iout,*) &
+ "Analytical (upper) and numerical (lower) gradient of theta"
+ write (iout,*)
+ do i=3,nres
+ do j=1,3
+ dcji=dc(j,i-2)
+ dc(j,i-2)=dcji+aincr
+ call chainbuild_cart
+ call int_from_cart1(.false.)
+ dthetanum(j,1,i)=(theta(i)-theta_s(i))/aincr
+ dc(j,i-2)=dcji
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dc(j,i-1)+aincr
+ call chainbuild_cart
+ dthetanum(j,2,i)=(theta(i)-theta_s(i))/aincr
+ dc(j,i-1)=dcji
+ enddo
+!el write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dtheta(j,1,i),j=1,3),&
+!el (dtheta(j,2,i),j=1,3)
+!el write (iout,'(5x,3f10.5,5x,3f10.5)') (dthetanum(j,1,i),j=1,3),&
+!el (dthetanum(j,2,i),j=1,3)
+!el write (iout,'(5x,3f10.5,5x,3f10.5)') &
+!el (dthetanum(j,1,i)/dtheta(j,1,i),j=1,3),&
+!el (dthetanum(j,2,i)/dtheta(j,2,i),j=1,3)
+!el write (iout,*)
+ enddo
+! Check gamma gradient
+ write (iout,*) &
+ "Analytical (upper) and numerical (lower) gradient of gamma"
+ do i=4,nres
+ do j=1,3
+ dcji=dc(j,i-3)
+ dc(j,i-3)=dcji+aincr
+ call chainbuild_cart
+ dphinum(j,1,i)=(phi(i)-phi_s(i))/aincr
+ dc(j,i-3)=dcji
+ dcji=dc(j,i-2)
+ dc(j,i-2)=dcji+aincr
+ call chainbuild_cart
+ dphinum(j,2,i)=(phi(i)-phi_s(i))/aincr
+ dc(j,i-2)=dcji
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dc(j,i-1)+aincr
+ call chainbuild_cart
+ dphinum(j,3,i)=(phi(i)-phi_s(i))/aincr
+ dc(j,i-1)=dcji
+ enddo
+!el write (iout,'(i5,3(3f10.5,5x))') i,(dphi(j,1,i),j=1,3),&
+!el (dphi(j,2,i),j=1,3),(dphi(j,3,i),j=1,3)
+!el write (iout,'(5x,3(3f10.5,5x))') (dphinum(j,1,i),j=1,3),&
+!el (dphinum(j,2,i),j=1,3),(dphinum(j,3,i),j=1,3)
+!el write (iout,'(5x,3(3f10.5,5x))') &
+!el (dphinum(j,1,i)/dphi(j,1,i),j=1,3),&
+!el (dphinum(j,2,i)/dphi(j,2,i),j=1,3),&
+!el (dphinum(j,3,i)/dphi(j,3,i),j=1,3)
+!el write (iout,*)
+ enddo
+! Check alpha gradient
+ write (iout,*) &
+ "Analytical (upper) and numerical (lower) gradient of alpha"
+ do i=2,nres-1
+ if(itype(i,1).ne.10) then
+ do j=1,3
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dcji+aincr
+ call chainbuild_cart
+ dalphanum(j,1,i)=(alph(i)-alph_s(i)) &
+ /aincr
+ dc(j,i-1)=dcji
+ dcji=dc(j,i)
+ dc(j,i)=dcji+aincr
+ call chainbuild_cart
+ dalphanum(j,2,i)=(alph(i)-alph_s(i)) &
+ /aincr
+ dc(j,i)=dcji
+ dcji=dc(j,i+nres)
+ dc(j,i+nres)=dc(j,i+nres)+aincr
+ call chainbuild_cart
+ dalphanum(j,3,i)=(alph(i)-alph_s(i)) &
+ /aincr
+ dc(j,i+nres)=dcji
+ enddo
+ endif
+!el write (iout,'(i5,3(3f10.5,5x))') i,(dalpha(j,1,i),j=1,3),&
+!el (dalpha(j,2,i),j=1,3),(dalpha(j,3,i),j=1,3)
+!el write (iout,'(5x,3(3f10.5,5x))') (dalphanum(j,1,i),j=1,3),&
+!el (dalphanum(j,2,i),j=1,3),(dalphanum(j,3,i),j=1,3)
+!el write (iout,'(5x,3(3f10.5,5x))') &
+!el (dalphanum(j,1,i)/dalpha(j,1,i),j=1,3),&
+!el (dalphanum(j,2,i)/dalpha(j,2,i),j=1,3),&
+!el (dalphanum(j,3,i)/dalpha(j,3,i),j=1,3)
+!el write (iout,*)
+ enddo
+! Check omega gradient
+ write (iout,*) &
+ "Analytical (upper) and numerical (lower) gradient of omega"
+ do i=2,nres-1
+ if(itype(i,1).ne.10) then
+ do j=1,3
+ dcji=dc(j,i-1)
+ dc(j,i-1)=dcji+aincr
+ call chainbuild_cart
+ domeganum(j,1,i)=(omeg(i)-omeg_s(i)) &
+ /aincr
+ dc(j,i-1)=dcji
+ dcji=dc(j,i)
+ dc(j,i)=dcji+aincr
+ call chainbuild_cart
+ domeganum(j,2,i)=(omeg(i)-omeg_s(i)) &
+ /aincr
+ dc(j,i)=dcji
+ dcji=dc(j,i+nres)
+ dc(j,i+nres)=dc(j,i+nres)+aincr
+ call chainbuild_cart
+ domeganum(j,3,i)=(omeg(i)-omeg_s(i)) &
+ /aincr
+ dc(j,i+nres)=dcji
+ enddo
+ endif
+!el write (iout,'(i5,3(3f10.5,5x))') i,(domega(j,1,i),j=1,3),&
+!el (domega(j,2,i),j=1,3),(domega(j,3,i),j=1,3)
+!el write (iout,'(5x,3(3f10.5,5x))') (domeganum(j,1,i),j=1,3),&
+!el (domeganum(j,2,i),j=1,3),(domeganum(j,3,i),j=1,3)
+!el write (iout,'(5x,3(3f10.5,5x))') &
+!el (domeganum(j,1,i)/domega(j,1,i),j=1,3),&
+!el (domeganum(j,2,i)/domega(j,2,i),j=1,3),&
+!el (domeganum(j,3,i)/domega(j,3,i),j=1,3)
+!el write (iout,*)
+ enddo
+ return
+ end subroutine checkintcartgrad
+!-----------------------------------------------------------------------------
+! q_measure.F
+!-----------------------------------------------------------------------------
+ real(kind=8) function qwolynes(seg1,seg2,flag,seg3,seg4)
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+ integer :: i,j,jl,k,l,il,kl,nl,np,ip,kp,seg1,seg2,seg3,seg4,secseg
+ integer :: kkk,nsep=3
+ real(kind=8) :: qm !dist,
+ real(kind=8) :: qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM,qqmax
+ logical :: lprn=.false.
+ logical :: flag
+! real(kind=8) :: sigm,x
+
+!el sigm(x)=0.25d0*x ! local function
+ qqmax=1.0d10
+ do kkk=1,nperm
+ qq = 0.0d0
+ nl=0
+ if(flag) then
+ do il=seg1+nsep,seg2
+ do jl=seg1,il-nsep
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2 + &
+ (cref(2,jl,kkk)-cref(2,il,kkk))**2 + &
+ (cref(3,jl,kkk)-cref(3,il,kkk))**2)
+ dij=dist(il,jl)
+ qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
+ nl=nl+1
+ d0ijCM=dsqrt( &
+ (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
+ (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
+ (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
+ dijCM=dist(il+nres,jl+nres)
+ qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
+ endif
+ qq = qq+qqij+qqijCM
+ enddo
+ enddo
+ qq = qq/nl
+ else
+ do il=seg1,seg2
+ if((seg3-il).lt.3) then
+ secseg=il+3
+ else
+ secseg=seg3
+ endif
+ do jl=secseg,seg4
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
+ (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
+ (cref(3,jl,kkk)-cref(3,il,kkk))**2)
+ dij=dist(il,jl)
+ qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
+ nl=nl+1
+ d0ijCM=dsqrt( &
+ (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
+ (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
+ (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
+ dijCM=dist(il+nres,jl+nres)
+ qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
+ endif
+ qq = qq+qqij+qqijCM
+ enddo
+ enddo
+ qq = qq/nl
+ endif
+ if (qqmax.le.qq) qqmax=qq
+ enddo
+ qwolynes=1.0d0-qqmax
+ return
+ end function qwolynes
+!-----------------------------------------------------------------------------
+ subroutine qwolynes_prim(seg1,seg2,flag,seg3,seg4)
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.MD'
+ integer :: i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,secseg
+ integer :: nsep=3, kkk
+!el real(kind=8) :: dist
+ real(kind=8) :: dij,d0ij,dijCM,d0ijCM
+ logical :: lprn=.false.
+ logical :: flag
+ real(kind=8) :: sim,dd0,fac,ddqij
+!el sigm(x)=0.25d0*x ! local function
+ do kkk=1,nperm
+ do i=0,nres
+ do j=1,3
+ dqwol(j,i)=0.0d0
+ dxqwol(j,i)=0.0d0
+ enddo
+ enddo
+ nl=0
+ if(flag) then
+ do il=seg1+nsep,seg2
+ do jl=seg1,il-nsep
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
+ (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
+ (cref(3,jl,kkk)-cref(3,il,kkk))**2)
+ dij=dist(il,jl)
+ sim = 1.0d0/sigm(d0ij)
+ sim = sim*sim
+ dd0 = dij-d0ij
+ fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il)-c(k,jl))*fac
+ dqwol(k,il)=dqwol(k,il)+ddqij
+ dqwol(k,jl)=dqwol(k,jl)-ddqij
+ enddo
+
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
+ nl=nl+1
+ d0ijCM=dsqrt( &
+ (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
+ (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
+ (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
+ dijCM=dist(il+nres,jl+nres)
+ sim = 1.0d0/sigm(d0ijCM)
+ sim = sim*sim
+ dd0=dijCM-d0ijCM
+ fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
+ dxqwol(k,il)=dxqwol(k,il)+ddqij
+ dxqwol(k,jl)=dxqwol(k,jl)-ddqij
+ enddo
+ endif
+ enddo
+ enddo
+ else
+ do il=seg1,seg2
+ if((seg3-il).lt.3) then
+ secseg=il+3
+ else
+ secseg=seg3
+ endif
+ do jl=secseg,seg4
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+ &
+ (cref(2,jl,kkk)-cref(2,il,kkk))**2+ &
+ (cref(3,jl,kkk)-cref(3,il,kkk))**2)
+ dij=dist(il,jl)
+ sim = 1.0d0/sigm(d0ij)
+ sim = sim*sim
+ dd0 = dij-d0ij
+ fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il)-c(k,jl))*fac
+ dqwol(k,il)=dqwol(k,il)+ddqij
+ dqwol(k,jl)=dqwol(k,jl)-ddqij
+ enddo
+ if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then
+ nl=nl+1
+ d0ijCM=dsqrt( &
+ (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ &
+ (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+ &
+ (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
+ dijCM=dist(il+nres,jl+nres)
+ sim = 1.0d0/sigm(d0ijCM)
+ sim=sim*sim
+ dd0 = dijCM-d0ijCM
+ fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
+ dxqwol(k,il)=dxqwol(k,il)+ddqij
+ dxqwol(k,jl)=dxqwol(k,jl)-ddqij
+ enddo
+ endif
+ enddo
+ enddo
+ endif
+ enddo
+ do i=0,nres
+ do j=1,3
+ dqwol(j,i)=dqwol(j,i)/nl
+ dxqwol(j,i)=dxqwol(j,i)/nl
+ enddo
+ enddo
+ return
+ end subroutine qwolynes_prim
+!-----------------------------------------------------------------------------
+ subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+ integer :: seg1,seg2,seg3,seg4
+ logical :: flag
+ real(kind=8),dimension(3,0:nres) :: qwolan,qwolxan
+ real(kind=8),dimension(3,0:2*nres) :: cdummy
+ real(kind=8) :: q1,q2
+ real(kind=8) :: delta=1.0d-10
+ integer :: i,j
+
+ do i=0,nres
+ do j=1,3
+ q1=qwolynes(seg1,seg2,flag,seg3,seg4)
+ cdummy(j,i)=c(j,i)
+ c(j,i)=c(j,i)+delta
+ q2=qwolynes(seg1,seg2,flag,seg3,seg4)
+ qwolan(j,i)=(q2-q1)/delta
+ c(j,i)=cdummy(j,i)
+ enddo
+ enddo
+ do i=0,nres
+ do j=1,3
+ q1=qwolynes(seg1,seg2,flag,seg3,seg4)
+ cdummy(j,i+nres)=c(j,i+nres)
+ c(j,i+nres)=c(j,i+nres)+delta
+ q2=qwolynes(seg1,seg2,flag,seg3,seg4)
+ qwolxan(j,i)=(q2-q1)/delta
+ c(j,i+nres)=cdummy(j,i+nres)
+ enddo
+ enddo
+! write(iout,*) "Numerical Q carteisan gradients backbone: "
+! do i=0,nct
+! write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
+! enddo
+! write(iout,*) "Numerical Q carteisan gradients side-chain: "
+! do i=0,nct
+! write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
+! enddo
+ return
+ end subroutine qwol_num
+!-----------------------------------------------------------------------------
+ subroutine EconstrQ
+! MD with umbrella_sampling using Wolyne's distance measure as a constraint
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.VAR'
+! include 'COMMON.MD'
+ use MD_data
+!#ifndef LANG0
+! include 'COMMON.LANGEVIN'
+!#else
+! include 'COMMON.LANGEVIN.lang0'
+!#endif
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+! include 'COMMON.TIME1'
+ real(kind=8) :: uzap1,uzap2,hm1,hm2,hmnum,ucdelan
+ real(kind=8),dimension(3,0:nres) :: dUcartan,dUxcartan,cdummy,&
+ duconst,duxconst
+ integer :: kstart,kend,lstart,lend,idummy
+ real(kind=8) :: delta=1.0d-7
+ integer :: i,j,k,ii
+ do i=0,nres
+ do j=1,3
+ duconst(j,i)=0.0d0
+ dudconst(j,i)=0.0d0
+ duxconst(j,i)=0.0d0
+ dudxconst(j,i)=0.0d0
+ enddo
+ enddo
+ Uconst=0.0d0
+ do i=1,nfrag
+ qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.,&
+ idummy,idummy)
+ Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
+! Calculating the derivatives of Constraint energy with respect to Q
+ Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),&
+ qinfrag(i,iset))
+! hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
+! hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
+! hmnum=(hm2-hm1)/delta
+! write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
+! & qinfrag(i,iset))
+! write(iout,*) "harmonicnum frag", hmnum
+! Calculating the derivatives of Q with respect to cartesian coordinates
+ call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.,&
+ idummy,idummy)
+! write(iout,*) "dqwol "
+! do ii=1,nres
+! write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
+! enddo
+! write(iout,*) "dxqwol "
+! do ii=1,nres
+! write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
+! enddo
+! Calculating numerical gradients of dU/dQi and dQi/dxi
+! call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
+! & ,idummy,idummy)
+! The gradients of Uconst in Cs
+ do ii=0,nres
+ do j=1,3
+ duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
+ dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
+ enddo
+ enddo
+ enddo
+ do i=1,npair
+ kstart=ifrag(1,ipair(1,i,iset),iset)
+ kend=ifrag(2,ipair(1,i,iset),iset)
+ lstart=ifrag(1,ipair(2,i,iset),iset)
+ lend=ifrag(2,ipair(2,i,iset),iset)
+ qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
+ Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
+! Calculating dU/dQ
+ Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
+! hm1=harmonic(qpair(i),qinpair(i,iset))
+! hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
+! hmnum=(hm2-hm1)/delta
+! write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
+! & qinpair(i,iset))
+! write(iout,*) "harmonicnum pair ", hmnum
+! Calculating dQ/dXi
+ call qwolynes_prim(kstart,kend,.false.,&
+ lstart,lend)
+! write(iout,*) "dqwol "
+! do ii=1,nres
+! write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
+! enddo
+! write(iout,*) "dxqwol "
+! do ii=1,nres
+! write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
+! enddo
+! Calculating numerical gradients
+! call qwol_num(kstart,kend,.false.
+! & ,lstart,lend)
+! The gradients of Uconst in Cs
+ do ii=0,nres
+ do j=1,3
+ duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
+ dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
+ enddo
+ enddo
+ enddo
+! write(iout,*) "Uconst inside subroutine ", Uconst
+! Transforming the gradients from Cs to dCs for the backbone
+ do i=0,nres
+ do j=i+1,nres
+ do k=1,3
+ dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
+ enddo
+ enddo
+ enddo
+! Transforming the gradients from Cs to dCs for the side chains
+ do i=1,nres
+ do j=1,3
+ dudxconst(j,i)=duxconst(j,i)
+ enddo
+ enddo
+! write(iout,*) "dU/ddc backbone "
+! do ii=0,nres
+! write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
+! enddo
+! write(iout,*) "dU/ddX side chain "
+! do ii=1,nres
+! write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
+! enddo
+! Calculating numerical gradients of dUconst/ddc and dUconst/ddx
+! call dEconstrQ_num
+ return
+ end subroutine EconstrQ
+!-----------------------------------------------------------------------------
+ subroutine dEconstrQ_num
+! Calculating numerical dUconst/ddc and dUconst/ddx
+! implicit real(kind=8) (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CONTROL'
+! include 'COMMON.VAR'
+! include 'COMMON.MD'
+ use MD_data
+!#ifndef LANG0
+! include 'COMMON.LANGEVIN'
+!#else
+! include 'COMMON.LANGEVIN.lang0'
+!#endif
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.NAMES'
+! include 'COMMON.TIME1'
+ real(kind=8) :: uzap1,uzap2
+ real(kind=8),dimension(3,0:nres) :: dUcartan,dUxcartan,cdummy
+ integer :: kstart,kend,lstart,lend,idummy
+ real(kind=8) :: delta=1.0d-7
+!el local variables
+ integer :: i,ii,j
+! real(kind=8) ::
+! For the backbone
+ do i=0,nres-1
+ do j=1,3
+ dUcartan(j,i)=0.0d0
+ cdummy(j,i)=dc(j,i)
+ dc(j,i)=dc(j,i)+delta
+ call chainbuild_cart
+ uzap2=0.0d0
+ do ii=1,nfrag
+ qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.,&
+ idummy,idummy)
+ uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),&
+ qinfrag(ii,iset))
+ enddo
+ do ii=1,npair
+ kstart=ifrag(1,ipair(1,ii,iset),iset)
+ kend=ifrag(2,ipair(1,ii,iset),iset)
+ lstart=ifrag(1,ipair(2,ii,iset),iset)
+ lend=ifrag(2,ipair(2,ii,iset),iset)
+ qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
+ uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),&
+ qinpair(ii,iset))
+ enddo
+ dc(j,i)=cdummy(j,i)
+ call chainbuild_cart
+ uzap1=0.0d0
+ do ii=1,nfrag
+ qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.,&
+ idummy,idummy)
+ uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),&
+ qinfrag(ii,iset))
+ enddo
+ do ii=1,npair
+ kstart=ifrag(1,ipair(1,ii,iset),iset)
+ kend=ifrag(2,ipair(1,ii,iset),iset)
+ lstart=ifrag(1,ipair(2,ii,iset),iset)
+ lend=ifrag(2,ipair(2,ii,iset),iset)
+ qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
+ uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),&
+ qinpair(ii,iset))
+ enddo
+ ducartan(j,i)=(uzap2-uzap1)/(delta)
+ enddo
+ enddo
+! Calculating numerical gradients for dU/ddx
+ do i=0,nres-1
+ duxcartan(j,i)=0.0d0
+ do j=1,3
+ cdummy(j,i)=dc(j,i+nres)
+ dc(j,i+nres)=dc(j,i+nres)+delta
+ call chainbuild_cart
+ uzap2=0.0d0
+ do ii=1,nfrag
+ qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.,&
+ idummy,idummy)
+ uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),&
+ qinfrag(ii,iset))
+ enddo
+ do ii=1,npair
+ kstart=ifrag(1,ipair(1,ii,iset),iset)
+ kend=ifrag(2,ipair(1,ii,iset),iset)
+ lstart=ifrag(1,ipair(2,ii,iset),iset)
+ lend=ifrag(2,ipair(2,ii,iset),iset)
+ qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
+ uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),&
+ qinpair(ii,iset))
+ enddo
+ dc(j,i+nres)=cdummy(j,i)
+ call chainbuild_cart
+ uzap1=0.0d0
+ do ii=1,nfrag
+ qfrag(ii)=qwolynes(ifrag(1,ii,iset),&
+ ifrag(2,ii,iset),.true.,idummy,idummy)
+ uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),&
+ qinfrag(ii,iset))
+ enddo
+ do ii=1,npair
+ kstart=ifrag(1,ipair(1,ii,iset),iset)
+ kend=ifrag(2,ipair(1,ii,iset),iset)
+ lstart=ifrag(1,ipair(2,ii,iset),iset)
+ lend=ifrag(2,ipair(2,ii,iset),iset)
+ qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
+ uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),&
+ qinpair(ii,iset))
+ enddo
+ duxcartan(j,i)=(uzap2-uzap1)/(delta)
+ enddo
+ enddo
+ write(iout,*) "Numerical dUconst/ddc backbone "
+ do ii=0,nres
+ write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
+ enddo
+! write(iout,*) "Numerical dUconst/ddx side-chain "
+! do ii=1,nres
+! write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
+! enddo
+ return
+ end subroutine dEconstrQ_num
+!-----------------------------------------------------------------------------
+! ssMD.F
+!-----------------------------------------------------------------------------
+ subroutine check_energies
+
+! use random, only: ran_number
+
+! implicit none
+! Includes
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.LOCAL'
+! include 'COMMON.GEO'
+
+! External functions
+!EL double precision ran_number
+!EL external ran_number
+
+! Local variables
+ integer :: i,j,k,l,lmax,p,pmax,countss
+ real(kind=8) :: rmin,rmax
+ real(kind=8) :: eij
+
+ real(kind=8) :: d
+ real(kind=8) :: wi,rij,tj,pj
+! return
+ countss=1
+ i=5
+ j=14
+
+ d=dsc(1)
+ rmin=2.0D0
+ rmax=12.0D0
+
+ lmax=10000
+ pmax=1
+
+ do k=1,3
+ c(k,i)=0.0D0
+ c(k,j)=0.0D0
+ c(k,nres+i)=0.0D0
+ c(k,nres+j)=0.0D0
+ enddo
+
+ do l=1,lmax
+
+!t wi=ran_number(0.0D0,pi)
+! wi=ran_number(0.0D0,pi/6.0D0)
+! wi=0.0D0
+!t tj=ran_number(0.0D0,pi)
+!t pj=ran_number(0.0D0,pi)
+! pj=ran_number(0.0D0,pi/6.0D0)
+! pj=0.0D0
+
+ do p=1,pmax
+!t rij=ran_number(rmin,rmax)
+
+ c(1,j)=d*sin(pj)*cos(tj)
+ c(2,j)=d*sin(pj)*sin(tj)
+ c(3,j)=d*cos(pj)
+
+ c(3,nres+i)=-rij
+
+ c(1,i)=d*sin(wi)
+ c(3,i)=-rij-d*cos(wi)
+
+ do k=1,3
+ dc(k,nres+i)=c(k,nres+i)-c(k,i)
+ dc_norm(k,nres+i)=dc(k,nres+i)/d
+ dc(k,nres+j)=c(k,nres+j)-c(k,j)
+ dc_norm(k,nres+j)=dc(k,nres+j)/d
+ enddo
+
+ call dyn_ssbond_ene(i,j,eij,countss)
+ enddo
+ enddo
+ call exit(1)
+ return
+ end subroutine check_energies
+!-----------------------------------------------------------------------------
+ subroutine dyn_ssbond_ene(resi,resj,eij,countss)
+! implicit none
+! Includes
+ use calc_data
+ use comm_sschecks
+! include 'DIMENSIONS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ use MD_data
+! include 'COMMON.MD'
+! use MD, only: totT,t_bath
+#endif
+#endif
+! External functions
+!EL double precision h_base
+!EL external h_base
+
+! Input arguments
+ integer :: resi,resj
+
+! Output arguments
+ real(kind=8) :: eij
+
+! Local variables
+ logical :: havebond
+ integer itypi,itypj,countss
+ real(kind=8) :: rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ real(kind=8) :: sig0ij,ljd,sig,fac,e1,e2
+ real(kind=8),dimension(3) :: dcosom1,dcosom2
+ real(kind=8) :: ed
+ real(kind=8) :: pom1,pom2
+ real(kind=8) :: ljA,ljB,ljXs
+ real(kind=8),dimension(1:3) :: d_ljB
+ real(kind=8) :: ssA,ssB,ssC,ssXs
+ real(kind=8) :: ssxm,ljxm,ssm,ljm
+ real(kind=8),dimension(1:3) :: d_ssxm,d_ljxm,d_ssm,d_ljm
+ real(kind=8) :: f1,f2,h1,h2,hd1,hd2
+ real(kind=8) :: omega,delta_inv,deltasq_inv,fac1,fac2
+!-------FIRST METHOD
+ real(kind=8) :: xm
+ real(kind=8),dimension(1:3) :: d_xm
+!-------END FIRST METHOD
+!-------SECOND METHOD
+!$$$ double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
+!-------END SECOND METHOD
+
+!-------TESTING CODE
+!el logical :: checkstop,transgrad
+!el common /sschecks/ checkstop,transgrad
+
+ integer :: icheck,nicheck,jcheck,njcheck
+ real(kind=8),dimension(-1:1) :: echeck
+ real(kind=8) :: deps,ssx0,ljx0
+!-------END TESTING CODE
+
+ eij=0.0d0
+ i=resi
+ j=resj
+
+!el allocate(dyn_ssbond_ij(iatsc_s:iatsc_e,nres))
+!el allocate(dyn_ssbond_ij(0:nres+4,nres))
+
+ itypi=itype(i,1)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+
+ itypj=itype(j,1)
+ xj=c(1,nres+j)-c(1,nres+i)
+ yj=c(2,nres+j)-c(2,nres+i)
+ zj=c(3,nres+j)-c(3,nres+i)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
+
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+! The following are set in sc_angular
+! erij(1)=xj*rij
+! erij(2)=yj*rij
+! erij(3)=zj*rij
+! om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
+! om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
+! om12=dxi*dxj+dyi*dyj+dzi*dzj
+ call sc_angular
+ rij=1.0D0/rij ! Reset this so it makes sense
+
+ sig0ij=sigma(itypi,itypj)
+ sig=sig0ij*dsqrt(1.0D0/sigsq)
+
+ ljXs=sig-sig0ij
+ ljA=eps1*eps2rt**2*eps3rt**2
+ ljB=ljA*bb_aq(itypi,itypj)
+ ljA=ljA*aa_aq(itypi,itypj)
+ ljxm=ljXs+(-2.0D0*aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0)
+
+ ssXs=d0cm
+ deltat1=1.0d0-om1
+ deltat2=1.0d0+om2
+ deltat12=om2-om1+2.0d0
+ cosphi=om12-om1*om2
+ ssA=akcm
+ ssB=akct*deltat12
+ ssC=ss_depth &
+ +akth*(deltat1*deltat1+deltat2*deltat2) &
+ +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
+ ssxm=ssXs-0.5D0*ssB/ssA
+
+!-------TESTING CODE
+!$$$c Some extra output
+!$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
+!$$$ ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+!$$$ ssx0=ssB*ssB-4.0d0*ssA*ssC
+!$$$ if (ssx0.gt.0.0d0) then
+!$$$ ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
+!$$$ else
+!$$$ ssx0=ssxm
+!$$$ endif
+!$$$ ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+!$$$ write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
+!$$$ & ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
+!$$$ return
+!-------END TESTING CODE
+
+!-------TESTING CODE
+! Stop and plot energy and derivative as a function of distance
+ if (checkstop) then
+ ssm=ssC-0.25D0*ssB*ssB/ssA
+ ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
+ if (ssm.lt.ljm .and. &
+ dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
+ nicheck=1000
+ njcheck=1
+ deps=0.5d-7
+ else
+ checkstop=.false.
+ endif
+ endif
+ if (.not.checkstop) then
+ nicheck=0
+ njcheck=-1
+ endif
+
+ do icheck=0,nicheck
+ do jcheck=-1,njcheck
+ if (checkstop) rij=(ssxm-1.0d0)+ &
+ ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
+!-------END TESTING CODE
+
+ if (rij.gt.ljxm) then
+ havebond=.false.
+ ljd=rij-ljXs
+ fac=(1.0D0/ljd)**expon
+ e1=fac*fac*aa_aq(itypi,itypj)
+ e2=fac*bb_aq(itypi,itypj)
+ eij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=eij*eps3rt
+ eps3der=eij*eps2rt
+ eij=eij*eps2rt*eps3rt
+
+ sigder=-sig/sigsq
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ ed=-expon*(e1+eij)/ljd
+ sigder=ed*sigder
+ eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
+ eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
+ eom12=eij*eps1_om12+eps2der*eps2rt_om12 &
+ -2.0D0*alf12*eps3der+sigder*sigsq_om12
+ else if (rij.lt.ssxm) then
+ havebond=.true.
+ ssd=rij-ssXs
+ eij=ssA*ssd*ssd+ssB*ssd+ssC
+
+ ed=2*akcm*ssd+akct*deltat12
+ pom1=akct*ssd
+ pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
+ eom1=-2*akth*deltat1-pom1-om2*pom2
+ eom2= 2*akth*deltat2+pom1-om1*pom2
+ eom12=pom2
+ else
+ omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi
+
+ d_ssxm(1)=0.5D0*akct/ssA
+ d_ssxm(2)=-d_ssxm(1)
+ d_ssxm(3)=0.0D0
+
+ d_ljxm(1)=sig0ij/sqrt(sigsq**3)
+ d_ljxm(2)=d_ljxm(1)*sigsq_om2
+ d_ljxm(3)=d_ljxm(1)*sigsq_om12
+ d_ljxm(1)=d_ljxm(1)*sigsq_om1
+
+!-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
+ xm=0.5d0*(ssxm+ljxm)
+ do k=1,3
+ d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
+ enddo
+ if (rij.lt.xm) then
+ havebond=.true.
+ ssm=ssC-0.25D0*ssB*ssB/ssA
+ d_ssm(1)=0.5D0*akct*ssB/ssA
+ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
+ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
+ d_ssm(3)=omega
+ f1=(rij-xm)/(ssxm-xm)
+ f2=(rij-ssxm)/(xm-ssxm)
+ h1=h_base(f1,hd1)
+ h2=h_base(f2,hd2)
+ eij=ssm*h1+Ht*h2
+ delta_inv=1.0d0/(xm-ssxm)
+ deltasq_inv=delta_inv*delta_inv
+ fac=ssm*hd1-Ht*hd2
+ fac1=deltasq_inv*fac*(xm-rij)
+ fac2=deltasq_inv*fac*(rij-ssxm)
+ ed=delta_inv*(Ht*hd2-ssm*hd1)
+ eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
+ eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
+ eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
+ else
+ havebond=.false.
+ ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
+ d_ljm(1)=-0.5D0*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)*ljB
+ d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
+ d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- &
+ alf12/eps3rt)
+ d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
+ f1=(rij-ljxm)/(xm-ljxm)
+ f2=(rij-xm)/(ljxm-xm)
+ h1=h_base(f1,hd1)
+ h2=h_base(f2,hd2)
+ eij=Ht*h1+ljm*h2
+ delta_inv=1.0d0/(ljxm-xm)
+ deltasq_inv=delta_inv*delta_inv
+ fac=Ht*hd1-ljm*hd2
+ fac1=deltasq_inv*fac*(ljxm-rij)
+ fac2=deltasq_inv*fac*(rij-xm)
+ ed=delta_inv*(ljm*hd2-Ht*hd1)
+ eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
+ eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
+ eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
+ endif
+!-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
+
+!-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
+!$$$ ssd=rij-ssXs
+!$$$ ljd=rij-ljXs
+!$$$ fac1=rij-ljxm
+!$$$ fac2=rij-ssxm
+!$$$
+!$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
+!$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
+!$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
+!$$$
+!$$$ ssm=ssC-0.25D0*ssB*ssB/ssA
+!$$$ d_ssm(1)=0.5D0*akct*ssB/ssA
+!$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
+!$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
+!$$$ d_ssm(3)=omega
+!$$$
+!$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
+!$$$ do k=1,3
+!$$$ d_ljm(k)=ljm*d_ljB(k)
+!$$$ enddo
+!$$$ ljm=ljm*ljB
+!$$$
+!$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC
+!$$$ d_ss(0)=2.0d0*ssA*ssd+ssB
+!$$$ d_ss(2)=akct*ssd
+!$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
+!$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
+!$$$ d_ss(3)=omega
+!$$$
+!$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj)
+!$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
+!$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1
+!$$$ do k=1,3
+!$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
+!$$$ & 2.0d0*ljB*fac1*d_ljxm(k))
+!$$$ enddo
+!$$$ ljf=ljm+ljf*ljB*fac1*fac1
+!$$$
+!$$$ f1=(rij-ljxm)/(ssxm-ljxm)
+!$$$ f2=(rij-ssxm)/(ljxm-ssxm)
+!$$$ h1=h_base(f1,hd1)
+!$$$ h2=h_base(f2,hd2)
+!$$$ eij=ss*h1+ljf*h2
+!$$$ delta_inv=1.0d0/(ljxm-ssxm)
+!$$$ deltasq_inv=delta_inv*delta_inv
+!$$$ fac=ljf*hd2-ss*hd1
+!$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
+!$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
+!$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
+!$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
+!$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
+!$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
+!$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
+!$$$
+!$$$ havebond=.false.
+!$$$ if (ed.gt.0.0d0) havebond=.true.
+!-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
+
+ endif
+
+ if (havebond) then
+!#ifndef CLUST
+!#ifndef WHAM
+! if (dyn_ssbond_ij(i,j).eq.1.0d300) then
+! write(iout,'(a15,f12.2,f8.1,2i5)')
+! & "SSBOND_E_FORM",totT,t_bath,i,j
+! endif
+!#endif
+!#endif
+ dyn_ssbond_ij(countss)=eij
+ else if (.not.havebond .and. dyn_ssbond_ij(countss).lt.1.0d300) then
+ dyn_ssbond_ij(countss)=1.0d300
+!#ifndef CLUST
+!#ifndef WHAM
+! write(iout,'(a15,f12.2,f8.1,2i5)')
+! & "SSBOND_E_BREAK",totT,t_bath,i,j
+!#endif
+!#endif
+ endif
+
+!-------TESTING CODE
+!el if (checkstop) then
+ if (jcheck.eq.0) write(iout,'(a,3f15.8,$)') &
+ "CHECKSTOP",rij,eij,ed
+ echeck(jcheck)=eij
+!el endif
+ enddo
+ if (checkstop) then
+ write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
+ endif
+ enddo
+ if (checkstop) then
+ transgrad=.true.
+ checkstop=.false.
+ endif
+!-------END TESTING CODE
+
+ do k=1,3
+ dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
+ dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
+ enddo
+ do k=1,3
+ gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ enddo
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k) &
+ +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ gvdwx(k,j)=gvdwx(k,j)+gg(k) &
+ +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ enddo
+!grad do k=i,j-1
+!grad do l=1,3
+!grad gvdwc(l,k)=gvdwc(l,k)+gg(l)
+!grad enddo
+!grad enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+
+ return
+ end subroutine dyn_ssbond_ene
+!--------------------------------------------------------------------------
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
+! implicit none
+! Includes
+ use calc_data
+ use comm_sschecks
+! include 'DIMENSIONS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ use MD_data
+! include 'COMMON.MD'
+! use MD, only: totT,t_bath
+#endif
+#endif
+ double precision h_base
+ external h_base
+
+!c Input arguments
+ integer resi,resj,resk,m,itypi,itypj,itypk
+
+!c Output arguments
+ double precision eij,eij1,eij2,eij3
+
+!c Local variables
+ logical havebond
+!c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+ double precision xik,yik,zik,xjk,yjk,zjk,dxk,dyk,dzk
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ eij=0.0
+ if (dtriss.eq.0) return
+ i=resi
+ j=resj
+ k=resk
+!C write(iout,*) resi,resj,resk
+ itypi=itype(i,1)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ itypj=itype(j,1)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+ itypk=itype(k,1)
+ xk=c(1,nres+k)
+ yk=c(2,nres+k)
+ zk=c(3,nres+k)
+ call to_box(xk,yk,zk)
+ dxk=dc_norm(1,nres+k)
+ dyk=dc_norm(2,nres+k)
+ dzk=dc_norm(3,nres+k)
+ dscj_inv=vbld_inv(k+nres)
+ xij=xj-xi
+ xik=xk-xi
+ xjk=xk-xj
+ yij=yj-yi
+ yik=yk-yi
+ yjk=yk-yj
+ zij=zj-zi
+ zik=zk-zi
+ zjk=zk-zj
+ rrij=(xij*xij+yij*yij+zij*zij)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+ rrik=(xik*xik+yik*yik+zik*zik)
+ rik=dsqrt(rrik)
+ rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+ rjk=dsqrt(rrjk)
+!C there are three combination of distances for each trisulfide bonds
+!C The first case the ith atom is the center
+!C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+!C distance y is second distance the a,b,c,d are parameters derived for
+!C this problem d parameter was set as a penalty currenlty set to 1.
+ if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
+ eij1=0.0d0
+ else
+ eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
+ endif
+!C second case jth atom is center
+ if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
+ eij2=0.0d0
+ else
+ eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
+ endif
+!C the third case kth atom is the center
+ if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
+ eij3=0.0d0
+ else
+ eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
+ endif
+!C eij2=0.0
+!C eij3=0.0
+!C eij1=0.0
+ eij=eij1+eij2+eij3
+!C write(iout,*)i,j,k,eij
+!C The energy penalty calculated now time for the gradient part
+!C derivative over rij
+ fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
+ gg(1)=xij*fac/rij
+ gg(2)=yij*fac/rij
+ gg(3)=zij*fac/rij
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,j)=gvdwx(m,j)+gg(m)
+ enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+!C now derivative over rik
+ fac=-eij1**2/dtriss* &
+ (-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xik*fac/rik
+ gg(2)=yik*fac/rik
+ gg(3)=zik*fac/rik
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+!C now derivative over rjk
+ fac=-eij2**2/dtriss* &
+ (-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- &
+ eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xjk*fac/rjk
+ gg(2)=yjk*fac/rjk
+ gg(3)=zjk*fac/rjk
+ do m=1,3
+ gvdwx(m,j)=gvdwx(m,j)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+ return
+ end subroutine triple_ssbond_ene
+
+
+
+!-----------------------------------------------------------------------------
+ real(kind=8) function h_base(x,deriv)
+! A smooth function going 0->1 in range [0,1]
+! It should NOT be called outside range [0,1], it will not work there.
+ implicit none
+
+! Input arguments
+ real(kind=8) :: x
+
+! Output arguments
+ real(kind=8) :: deriv
+
+! Local variables
+ real(kind=8) :: xsq
+
+
+! Two parabolas put together. First derivative zero at extrema
+!$$$ if (x.lt.0.5D0) then
+!$$$ h_base=2.0D0*x*x
+!$$$ deriv=4.0D0*x
+!$$$ else
+!$$$ deriv=1.0D0-x
+!$$$ h_base=1.0D0-2.0D0*deriv*deriv
+!$$$ deriv=4.0D0*deriv
+!$$$ endif
+
+! Third degree polynomial. First derivative zero at extrema
+ h_base=x*x*(3.0d0-2.0d0*x)
+ deriv=6.0d0*x*(1.0d0-x)
+
+! Fifth degree polynomial. First and second derivatives zero at extrema
+!$$$ xsq=x*x
+!$$$ h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
+!$$$ deriv=x-1.0d0
+!$$$ deriv=deriv*deriv
+!$$$ deriv=30.0d0*xsq*deriv
+
+ return
+ end function h_base
+!-----------------------------------------------------------------------------
+ subroutine dyn_set_nss
+! Adjust nss and other relevant variables based on dyn_ssbond_ij
+! implicit none
+ use MD_data, only: totT,t_bath
+! Includes
+! include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.SETUP'
+! include 'COMMON.MD'
+! Local variables
+ real(kind=8) :: emin
+ integer :: i,j,imin,ierr,k
+ integer :: diff,allnss,newnss
+ integer,dimension(maxdim) :: allflag,allihpb,alljhpb,& !(maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
+ newihpb,newjhpb,aliass
+ logical :: found
+ integer,dimension(0:nfgtasks) :: i_newnss
+ integer,dimension(0:nfgtasks) :: displ
+ integer,dimension(maxdim) :: g_newihpb,g_newjhpb !(maxdim)(maxdim=(maxres-1)*(maxres-2)/2)
+ integer :: g_newnss
+
+ allnss=0
+ k=0
+ do i=1,nres-1
+ do j=i+1,nres
+ if ((itype(i,1).eq.1).and.(itype(j,1).eq.1)) then
+ k=k+1
+ if (dyn_ssbond_ij(k).lt.1.0d300) then
+ allnss=allnss+1
+ allflag(allnss)=0
+ allihpb(allnss)=i
+ alljhpb(allnss)=j
+ aliass(allnss)=k
+ endif
+ endif
+ enddo
+ enddo
+
+!mc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ 1 emin=1.0d300
+ do i=1,allnss
+ if (allflag(i).eq.0 .and. &
+ dyn_ssbond_ij(aliass(allnss)).lt.emin) then
+ emin=dyn_ssbond_ij(aliass(allnss))
+ imin=i
+ endif
+ enddo
+ if (emin.lt.1.0d300) then
+ allflag(imin)=1
+ do i=1,allnss
+ if (allflag(i).eq.0 .and. &
+ (allihpb(i).eq.allihpb(imin) .or. &
+ alljhpb(i).eq.allihpb(imin) .or. &
+ allihpb(i).eq.alljhpb(imin) .or. &
+ alljhpb(i).eq.alljhpb(imin))) then
+ allflag(i)=-1
+ endif
+ enddo
+ goto 1
+ endif
+
+!mc write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)
+
+ newnss=0
+ do i=1,allnss
+ if (allflag(i).eq.1) then
+ newnss=newnss+1
+ newihpb(newnss)=allihpb(i)
+ newjhpb(newnss)=alljhpb(i)
+ endif
+ enddo
+
+#ifdef MPI
+ if (nfgtasks.gt.1)then
+
+ call MPI_Reduce(newnss,g_newnss,1,&
+ MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
+ call MPI_Gather(newnss,1,MPI_INTEGER,&
+ i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
+ displ(0)=0
+ do i=1,nfgtasks-1,1
+ displ(i)=i_newnss(i-1)+displ(i-1)
+ enddo
+ call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,&
+ g_newihpb,i_newnss,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,&
+ g_newjhpb,i_newnss,displ,MPI_INTEGER,&
+ king,FG_COMM,IERR)
+ if(fg_rank.eq.0) then
+! print *,'g_newnss',g_newnss
+! print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
+! print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
+ newnss=g_newnss
+ do i=1,newnss
+ newihpb(i)=g_newihpb(i)
+ newjhpb(i)=g_newjhpb(i)
+ enddo
+ endif
+ endif
+#endif
+
+ diff=newnss-nss
+
+!mc write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
+! print *,newnss,nss,maxdim
+ do i=1,nss
+ found=.false.
+! print *,newnss
+ do j=1,newnss
+!! print *,j
+ if (idssb(i).eq.newihpb(j) .and. &
+ jdssb(i).eq.newjhpb(j)) found=.true.
+ enddo
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
+! write(iout,*) "found",found,i,j
+ if (.not.found.and.fg_rank.eq.0) &
+ write(iout,'(a15,f12.2,f8.1,2i5)') &
+ "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
+#endif
+ enddo
+
+ do i=1,newnss
+ found=.false.
+ do j=1,nss
+! print *,i,j
+ if (newihpb(i).eq.idssb(j) .and. &
+ newjhpb(i).eq.jdssb(j)) found=.true.
+ enddo
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
+! write(iout,*) "found",found,i,j
+ if (.not.found.and.fg_rank.eq.0) &
+ write(iout,'(a15,f12.2,f8.1,2i5)') &
+ "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
+#endif
+ enddo
+!#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
+ nss=newnss
+ do i=1,nss
+ idssb(i)=newihpb(i)
+ jdssb(i)=newjhpb(i)
+ enddo
+!#else
+! nss=0
+!#endif
+
+ return
+ end subroutine dyn_set_nss
+! Lipid transfer energy function
+ subroutine Eliptransfer(eliptran)
+!C this is done by Adasko
+!C print *,"wchodze"
+!C structure of box:
+!C water
+!C--bordliptop-- buffore starts
+!C--bufliptop--- here true lipid starts
+!C lipid
+!C--buflipbot--- lipid ends buffore starts
+!C--bordlipbot--buffore ends
+ real(kind=8) :: fracinbuf,eliptran,sslip,positi,ssgradlip
+ integer :: i
+ eliptran=0.0
+! print *, "I am in eliptran"
+ do i=ilip_start,ilip_end
+!C do i=1,1
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1).or.(i.eq.nres))&
+ cycle
+
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
+ if (positi.le.0.0) positi=positi+boxzsize
+!C print *,i
+!C first for peptide groups
+!c for each residue check if it is in lipid or lipid water border area
+ if ((positi.gt.bordlipbot) &
+ .and.(positi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (positi.lt.buflipbot) then
+!C what fraction I am in
+ fracinbuf=1.0d0- &
+ ((positi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+!C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+
+!C print *,"doing sccale for lower part"
+!C print *,i,sslip,fracinbuf,ssgradlip
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+!C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+!C print *, "doing sscalefor top part"
+!C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+!C print *,"I am in true lipid"
+ endif
+!C else
+!C eliptran=elpitran+0.0 ! I am in water
+ endif
+ if (energy_dec) write(iout,*) i,"eliptran=",eliptran,positi,sslip
+ enddo
+! here starts the side chain transfer
+ do i=ilip_start,ilip_end
+ if (itype(i,1).eq.ntyp1) cycle
+ positi=(mod(c(3,i+nres),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+!C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
+!c for each residue check if it is in lipid or lipid water border area
+!C respos=mod(c(3,i+nres),boxzsize)
+!C print *,positi,bordlipbot,buflipbot
+ if ((positi.gt.bordlipbot) &
+ .and.(positi.lt.bordliptop)) then
+!C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0- &
+ ((positi-bordlipbot)/lipbufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
+ gliptranx(3,i)=gliptranx(3,i) &
+ +ssgradlip*liptranene(itype(i,1))
+ gliptranc(3,i-1)= gliptranc(3,i-1) &
+ +ssgradlip*liptranene(itype(i,1))
+!C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0- &
+ ((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i,1))
+ gliptranx(3,i)=gliptranx(3,i) &
+ +ssgradlip*liptranene(itype(i,1))
+ gliptranc(3,i-1)= gliptranc(3,i-1) &
+ +ssgradlip*liptranene(itype(i,1))
+!C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i,1))
+!C print *,"I am in true lipid"
+ endif
+ endif ! if in lipid or buffor
+!C else
+!C eliptran=elpitran+0.0 ! I am in water
+ if (energy_dec) write(iout,*) i,"eliptran=",eliptran
+ enddo
+ return
+ end subroutine Eliptransfer
+!----------------------------------NANO FUNCTIONS
+!C-----------------------------------------------------------------------
+!C-----------------------------------------------------------
+!C This subroutine is to mimic the histone like structure but as well can be
+!C utilizet to nanostructures (infinit) small modification has to be used to
+!C make it finite (z gradient at the ends has to be changes as well as the x,y
+!C gradient has to be modified at the ends
+!C The energy function is Kihara potential
+!C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+!C 4eps is depth of well sigma is r_minimum r is distance from center of tube
+!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
+!C simple Kihara potential
+ subroutine calctube(Etube)
+ real(kind=8),dimension(3) :: vectube
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, &
+ sc_aa_tube,sc_bb_tube
+ integer :: i,j,iti
+ Etube=0.0d0
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+!C first we calculate the distance from tube center
+!C for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+ xmin=boxxsize
+ ymin=boxysize
+! Find minimum distance in periodic box
+ do j=-1,1
+ vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6- &
+ 6.0d0*pep_bb_tube)/rdiff6/rdiff
+!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+!C &rdiff,fac
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ enddo
+!C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+!C print *,gg_tube(1,0),"TU"
+
+
+ do i=itube_start,itube_end
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i,1)
+!C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1) &
+!C in UNRES uncomment the line below as GLY has no side-chain...
+!C .or.(iti.eq.10)
+ ) cycle
+ xmin=boxxsize
+ ymin=boxysize
+ do j=-1,1
+ vectube(1)=mod((c(1,i+nres)),boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i+nres)),boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+!C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2),
+!C & tubecenter(2)
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+ fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- &
+ 6.0d0*sc_bb_tube/rdiff6/rdiff
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ enddo
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)
+ enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calctube
+!C TO DO 1) add to total energy
+!C 2) add to gradient summation
+!C 3) add reading parameters (AND of course oppening of PARAM file)
+!C 4) add reading the center of tube
+!C 5) add COMMONs
+!C 6) add to zerograd
+!C 7) allocate matrices
+
+
+!C-----------------------------------------------------------------------
+!C-----------------------------------------------------------
+!C This subroutine is to mimic the histone like structure but as well can be
+!C utilizet to nanostructures (infinit) small modification has to be used to
+!C make it finite (z gradient at the ends has to be changes as well as the x,y
+!C gradient has to be modified at the ends
+!C The energy function is Kihara potential
+!C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6)
+!C 4eps is depth of well sigma is r_minimum r is distance from center of tube
+!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a
+!C simple Kihara potential
+ subroutine calctube2(Etube)
+ real(kind=8),dimension(3) :: vectube
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,&
+ sstube,ssgradtube,sc_aa_tube,sc_bb_tube
+ integer:: i,j,iti
+ Etube=0.0d0
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+!C first we calculate the distance from tube center
+!C first sugare-phosphate group for NARES this would be peptide group
+!C for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+!C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+!C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+!C vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+!C if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+ xmin=boxxsize
+ ymin=boxysize
+ do j=-1,1
+ vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
+ vectube(1)=vectube(1)+boxxsize*j
+ vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+ vectube(2)=vectube(2)+boxysize*j
+
+ xminact=abs(vectube(1)-tubecenter(1))
+ yminact=abs(vectube(2)-tubecenter(2))
+ if (xmin.gt.xminact) then
+ xmin=xminact
+ xtemp=vectube(1)
+ endif
+ if (ymin.gt.yminact) then
+ ymin=yminact
+ ytemp=vectube(2)
+ endif
+ enddo
+ vectube(1)=xtemp
+ vectube(2)=ytemp
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C THIS FRAGMENT MAKES TUBE FINITE
+ positi=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize)
+ if (positi.le.0) positi=positi+boxzsize
+!C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
+!c for each residue check if it is in lipid or lipid water border area
+!C respos=mod(c(3,i+nres),boxzsize)
+!C print *,positi,bordtubebot,buftubebot,bordtubetop
+ if ((positi.gt.bordtubebot) &
+ .and.(positi.lt.bordtubetop)) then
+!C the energy transfer exist
+ if (positi.lt.buftubebot) then
+ fracinbuf=1.0d0- &
+ ((positi-bordtubebot)/tubebufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C print *,"doing sccale for lower part"
+ elseif (positi.gt.buftubetop) then
+ fracinbuf=1.0d0- &
+ ((bordtubetop-positi)/tubebufthick)
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ sstube=1.0d0
+ ssgradtube=0.0d0
+ enetube(i)=enetube(i)+sstube*tubetranenepep
+!C print *,"I am in true lipid"
+ endif
+ else
+!C sstube=0.0d0
+!C ssgradtube=0.0d0
+ cycle
+ endif ! if in lipid or buffor
+
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=enetube(i)+sstube* &
+ (pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6)
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6- &
+ 6.0d0*pep_bb_tube)/rdiff6/rdiff*sstube
+!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+!C &rdiff,fac
+
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ gg_tube(3,i)=gg_tube(3,i) &
+ +ssgradtube*enetube(i)/sstube/2.0d0
+ gg_tube(3,i-1)= gg_tube(3,i-1) &
+ +ssgradtube*enetube(i)/sstube/2.0d0
+
+ enddo
+!C basically thats all code now we split for side-chains (REMEMBER to sum up at the END)
+!C print *,gg_tube(1,0),"TU"
+ do i=itube_start,itube_end
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i,1)
+!C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1) &
+!!C in UNRES uncomment the line below as GLY has no side-chain...
+ .or.(iti.eq.10) &
+ ) cycle
+ vectube(1)=c(1,i+nres)
+ vectube(1)=mod(vectube(1),boxxsize)
+ if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
+ vectube(2)=c(2,i+nres)
+ vectube(2)=mod(vectube(2),boxysize)
+ if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
+
+ vectube(1)=vectube(1)-tubecenter(1)
+ vectube(2)=vectube(2)-tubecenter(2)
+!C THIS FRAGMENT MAKES TUBE FINITE
+ positi=(mod(c(3,i+nres),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+!C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
+!c for each residue check if it is in lipid or lipid water border area
+!C respos=mod(c(3,i+nres),boxzsize)
+!C print *,positi,bordtubebot,buftubebot,bordtubetop
+
+ if ((positi.gt.bordtubebot) &
+ .and.(positi.lt.bordtubetop)) then
+!C the energy transfer exist
+ if (positi.lt.buftubebot) then
+ fracinbuf=1.0d0- &
+ ((positi-bordtubebot)/tubebufthick)
+!C lipbufthick is thickenes of lipid buffore
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
+!C print *,ssgradtube, sstube,tubetranene(itype(i,1))
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C print *,"doing sccale for lower part"
+ elseif (positi.gt.buftubetop) then
+ fracinbuf=1.0d0- &
+ ((bordtubetop-positi)/tubebufthick)
+
+ sstube=sscalelip(fracinbuf)
+ ssgradtube=sscagradlip(fracinbuf)/tubebufthick
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
+!C gg_tube_SC(3,i)=gg_tube_SC(3,i)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C gg_tube(3,i-1)= gg_tube(3,i-1)
+!C &+ssgradtube*tubetranene(itype(i,1))
+!C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ sstube=1.0d0
+ ssgradtube=0.0d0
+ enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1))
+!C print *,"I am in true lipid"
+ endif
+ else
+!C sstube=0.0d0
+!C ssgradtube=0.0d0
+ cycle
+ endif ! if in lipid or buffor
+!CEND OF FINITE FRAGMENT
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+ vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6)&
+ *sstube+enetube(i+nres)
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff-&
+ 6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube
+!C now direction of gg_tube vector
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ gg_tube_SC(3,i)=gg_tube_SC(3,i) &
+ +ssgradtube*enetube(i+nres)/sstube
+ gg_tube(3,i-1)= gg_tube(3,i-1) &
+ +ssgradtube*enetube(i+nres)/sstube
+
+ enddo
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)
+ enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calctube2
+!=====================================================================================================================================
+ subroutine calcnano(Etube)
+ use MD_data, only:totTafm
+ real(kind=8),dimension(3) :: vectube,cm
+
+ real(kind=8) :: Etube,xtemp,xminact,yminact,&
+ ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,&
+ sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact,tubezcenter,xi,yi,zi!,&
+! vecsim,vectrue
+ real(kind=8) :: eps,sig,aa_tub_lip,bb_tub_lip
+ integer:: i,j,iti,r,ilol,ityp
+! totTafm=2.0
+ Etube=0.0d0
+ call to_box(tubecenter(1),tubecenter(2),tubecenter(3))
+! print *,itube_start,itube_end,"poczatek"
+ do i=itube_start,itube_end
+ enetube(i)=0.0d0
+ enetube(i+nres)=0.0d0
+ enddo
+!C first we calculate the distance from tube center
+!C first sugare-phosphate group for NARES this would be peptide group
+!C for UNRES
+ do i=itube_start,itube_end
+!C lets ommit dummy atoms for now
+ if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle
+!C now calculate distance from center of tube and direction vectors
+
+! do j=-1,1
+ xi=(c(1,i)+c(1,i+1))/2.0d0
+ yi=(c(2,i)+c(2,i+1))/2.0d0
+ zi=((c(3,i)+c(3,i+1))/2.0d0)
+ call to_box(xi,yi,zi)
+! tubezcenter=totTafm*velNANOconst+tubecenter(3)
+
+ vectube(1)=boxshift(xi-tubecenter(1),boxxsize)
+ vectube(2)=boxshift(yi-tubecenter(2),boxysize)
+ vectube(3)=boxshift(zi-tubecenter(3),boxzsize)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+!C vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*pep_aa_tube/rdiff6- &
+ 6.0d0*pep_bb_tube)/rdiff6/rdiff
+!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i),
+!C &rdiff,fac
+ if (acavtubpep.eq.0.0d0) then
+!C go to 667
+ enecavtube(i)=0.0
+ faccav=0.0
+ else
+ denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6)
+ enecavtube(i)= &
+ (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep) &
+ /denominator
+ enecavtube(i)=0.0
+ faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff)) &
+ *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff) &
+ +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) &
+ /denominator**2.0d0
+!C faccav=0.0
+!C fac=fac+faccav
+!C 667 continue
+ endif
+ if (energy_dec) write(iout,*),"ETUBE_PEP",i,rdiff,enetube(i),enecavtube(i)
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0
+ gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0
+ enddo
+ enddo
+
+ do i=itube_start,itube_end
+ enecavtube(i)=0.0d0
+!C Lets not jump over memory as we use many times iti
+ iti=itype(i,1)
+!C lets ommit dummy atoms for now
+ if ((iti.eq.ntyp1) &
+!C in UNRES uncomment the line below as GLY has no side-chain...
+!C .or.(iti.eq.10)
+ ) cycle
+ xi=c(1,i+nres)
+ yi=c(2,i+nres)
+ zi=c(3,i+nres)
+ call to_box(xi,yi,zi)
+ tubezcenter=totTafm*velNANOconst+tubecenter(3)
+
+ vectube(1)=boxshift(xi-tubecenter(1),boxxsize)
+ vectube(2)=boxshift(yi-tubecenter(2),boxysize)
+ vectube(3)=boxshift(zi-tubecenter(3),boxzsize)
+
+
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+ sc_aa_tube=sc_aa_tube_par(iti)
+ sc_bb_tube=sc_bb_tube_par(iti)
+ enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+!C enetube(i+nres)=0.0d0
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- &
+ 6.0d0*sc_bb_tube/rdiff6/rdiff
+!C fac=0.0
+!C now direction of gg_tube vector
+!C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12)
+ if (acavtub(iti).eq.0.0d0) then
+!C go to 667
+ enecavtube(i+nres)=0.0d0
+ faccav=0.0d0
+ else
+ denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6)
+ enecavtube(i+nres)= &
+ (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti)) &
+ /denominator
+!C enecavtube(i)=0.0
+ faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff)) &
+ *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff) &
+ +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0) &
+ /denominator**2.0d0
+!C faccav=0.0
+ fac=fac+faccav
+!C 667 continue
+ endif
+!C print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator,
+!C & enecavtube(i),faccav
+!C print *,"licz=",
+!C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti))
+!C print *,"finene=",enetube(i+nres)+enecavtube(i)
+ do j=1,3
+ gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ if (energy_dec) write(iout,*),"ETUBE",i,rdiff,enetube(i+nres),enecavtube(i+nres)
+ enddo
+
+
+
+ do i=itube_start,itube_end
+ Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) &
+ +enecavtube(i+nres)
+ enddo
+
+ do i=ilipbond_start_tub,ilipbond_end_tub
+ ityp=itype(i,4)
+! print *,"ilipbond_start",ilipbond_start,i,ityp
+ if (ityp.gt.ntyp_molec(4)) cycle
+!C now calculate distance from center of tube and direction vectors
+ eps=lip_sig(ityp,18)*4.0d0
+ sig=lip_sig(ityp,18)
+ aa_tub_lip=eps/(sig**12)
+ bb_tub_lip=eps/(sig**6)
+! do j=-1,1
+ xi=c(1,i)
+ yi=c(2,i)
+ zi=c(3,i)
+ call to_box(xi,yi,zi)
+! tubezcenter=totTafm*velNANOconst+tubecenter(3)
+
+ vectube(1)=boxshift(xi-tubecenter(1),boxxsize)
+ vectube(2)=boxshift(yi-tubecenter(2),boxysize)
+ vectube(3)=boxshift(zi-tubecenter(3),boxzsize)
+
+!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1)
+!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2)
+!C as the tube is infinity we do not calculate the Z-vector use of Z
+!C as chosen axis
+!C vectube(3)=0.0d0
+!C now calculte the distance
+ tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2)
+!C now normalize vector
+ vectube(1)=vectube(1)/tub_r
+ vectube(2)=vectube(2)/tub_r
+ vectube(3)=vectube(3)/tub_r
+!C calculte rdiffrence between r and r0
+ rdiff=tub_r-tubeR0
+!C and its 6 power
+ rdiff6=rdiff**6.0d0
+!C for vectorization reasons we will sumup at the end to avoid depenence of previous
+ enetube(i)=aa_tub_lip/rdiff6**2.0d0+bb_tub_lip/rdiff6
+ Etube=Etube+enetube(i)
+!C write(iout,*) "TU13",i,rdiff6,enetube(i)
+!C print *,rdiff,rdiff6,pep_aa_tube
+!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6
+!C now we calculate gradient
+ fac=(-12.0d0*aa_tub_lip/rdiff6- &
+ 6.0d0*bb_tub_lip)/rdiff6/rdiff
+ do j=1,3
+ gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac
+ enddo
+ if (energy_dec) write(iout,*) "ETUBLIP",i,rdiff,enetube(i+nres)
+ enddo
+
+
+!-----------------------------------------------------------------------
+ if (fg_rank.eq.0) then
+ if (velNANOconst.ne.0) then
+ do j=1,3
+ cm(j)=0.0d0
+ enddo
+ do i=1,inanomove
+ ilol=inanotab(i)
+ do j=1,3
+ cm(j)=cm(j)+c(j,ilol)
+ enddo
+ enddo
+ do j=1,3
+ cm(j)=cm(j)/inanomove
+ enddo
+ vecsim=velNANOconst*totTafm+distnanoinit
+ vectrue=cm(3)-tubecenter(3)
+ etube=etube+0.5d0*forcenanoconst*( vectrue-vecsim)**2
+ fac=forcenanoconst*(vectrue-vecsim)/inanomove
+ do i=1,inanomove
+ ilol=inanotab(i)
+ gg_tube(3,ilol-1)=gg_tube(3,ilol-1)+fac
+ enddo
+ endif
+ endif
+! do i=1,20
+! print *,"begin", i,"a"
+! do r=1,10000
+! rdiff=r/100.0d0
+! rdiff6=rdiff**6.0d0
+! sc_aa_tube=sc_aa_tube_par(i)
+! sc_bb_tube=sc_bb_tube_par(i)
+! enetube(i)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6
+! denominator=(1.0d0+dcavtub(i)*rdiff6*rdiff6)
+! enecavtube(i)= &
+! (bcavtub(i)*rdiff+acavtub(i)*dsqrt(rdiff)+ccavtub(i)) &
+! /denominator
+
+! print '(5(f10.3,1x))',rdiff,enetube(i),enecavtube(i),enecavtube(i)+enetube(i)
+! enddo
+! print *,"end",i,"a"
+! enddo
+!C print *,"ETUBE", etube
+ return
+ end subroutine calcnano
+
+!===============================================
+!--------------------------------------------------------------------------------
+!C first for shielding is setting of function of side-chains
+
+ subroutine set_shield_fac2
+ real(kind=8) :: div77_81=0.974996043d0, &
+ div4_81=0.2222222222d0
+ real (kind=8) :: dist_pep_side,dist_side_calf,dist_pept_group, &
+ scale_fac_dist,fac_help_scale,VofOverlap,VolumeTotal,costhet,&
+ short,long,sinthet,costhet_fac,sh_frac_dist,rkprim,cosphi, &
+ sinphi,cosphi_fac,pep_side0pept_group,cosalfa,fac_alfa_sin
+!C the vector between center of side_chain and peptide group
+ real(kind=8),dimension(3) :: pep_side_long,side_calf, &
+ pept_group,costhet_grad,cosphi_grad_long, &
+ cosphi_grad_loc,pep_side_norm,side_calf_norm, &
+ sh_frac_dist_grad,pep_side
+ integer i,j,k
+!C write(2,*) "ivec",ivec_start,ivec_end
+ do i=1,nres
+ fac_shield(i)=0.0d0
+ ishield_list(i)=0
+ do j=1,3
+ grad_shield(j,i)=0.0d0
+ enddo
+ enddo
+ do i=ivec_start,ivec_end
+!C do i=1,nres-1
+!C if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
+! ishield_list(i)=0
+ if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle
+!Cif there two consequtive dummy atoms there is no peptide group between them
+!C the line below has to be changed for FGPROC>1
+ VolumeTotal=0.0
+ do k=1,nres
+ if ((itype(k,1).eq.ntyp1).or.(itype(k,1).eq.10)) cycle
+ dist_pep_side=0.0
+ dist_side_calf=0.0
+ do j=1,3
+!C first lets set vector conecting the ithe side-chain with kth side-chain
+ pep_side(j)=c(j,k+nres)-(c(j,i)+c(j,i+1))/2.0d0
+!C pep_side(j)=2.0d0
+!C and vector conecting the side-chain with its proper calfa
+ side_calf(j)=c(j,k+nres)-c(j,k)
+!C side_calf(j)=2.0d0
+ pept_group(j)=c(j,i)-c(j,i+1)
+!C lets have their lenght
+ dist_pep_side=pep_side(j)**2+dist_pep_side
+ dist_side_calf=dist_side_calf+side_calf(j)**2
+ dist_pept_group=dist_pept_group+pept_group(j)**2
+ enddo
+ dist_pep_side=sqrt(dist_pep_side)
+ dist_pept_group=sqrt(dist_pept_group)
+ dist_side_calf=sqrt(dist_side_calf)
+ do j=1,3
+ pep_side_norm(j)=pep_side(j)/dist_pep_side
+ side_calf_norm(j)=dist_side_calf
+ enddo
+!C now sscale fraction
+ sh_frac_dist=-(dist_pep_side-rpp(1,1)-buff_shield)/buff_shield
+! print *,buff_shield,"buff",sh_frac_dist
+!C now sscale
+ if (sh_frac_dist.le.0.0) cycle
+!C print *,ishield_list(i),i
+!C If we reach here it means that this side chain reaches the shielding sphere
+!C Lets add him to the list for gradient
+ ishield_list(i)=ishield_list(i)+1
+!C ishield_list is a list of non 0 side-chain that contribute to factor gradient
+!C this list is essential otherwise problem would be O3
+ shield_list(ishield_list(i),i)=k
+!C Lets have the sscale value
+ if (sh_frac_dist.gt.1.0) then
+ scale_fac_dist=1.0d0
+ do j=1,3
+ sh_frac_dist_grad(j)=0.0d0
+ enddo
+ else
+ scale_fac_dist=-sh_frac_dist*sh_frac_dist &
+ *(2.0d0*sh_frac_dist-3.0d0)
+ fac_help_scale=6.0d0*(sh_frac_dist-sh_frac_dist**2) &
+ /dist_pep_side/buff_shield*0.5d0
+ do j=1,3
+ sh_frac_dist_grad(j)=fac_help_scale*pep_side(j)
+!C sh_frac_dist_grad(j)=0.0d0
+!C scale_fac_dist=1.0d0
+!C print *,"jestem",scale_fac_dist,fac_help_scale,
+!C & sh_frac_dist_grad(j)
+ enddo
+ endif
+!C this is what is now we have the distance scaling now volume...
+ short=short_r_sidechain(itype(k,1))
+ long=long_r_sidechain(itype(k,1))
+ costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2)
+ sinthet=short/dist_pep_side*costhet
+! print *,"SORT",short,long,sinthet,costhet
+!C now costhet_grad
+!C costhet=0.6d0
+!C sinthet=0.8
+ costhet_fac=costhet**3*short**2*(-0.5d0)/dist_pep_side**4
+!C sinthet_fac=costhet**2*0.5d0*(short**3/dist_pep_side**4*costhet
+!C & -short/dist_pep_side**2/costhet)
+!C costhet_fac=0.0d0
+ do j=1,3
+ costhet_grad(j)=costhet_fac*pep_side(j)
+ enddo
+!C remember for the final gradient multiply costhet_grad(j)
+!C for side_chain by factor -2 !
+!C fac alfa is angle between CB_k,CA_k, CA_i,CA_i+1
+!C pep_side0pept_group is vector multiplication
+ pep_side0pept_group=0.0d0
+ do j=1,3
+ pep_side0pept_group=pep_side0pept_group+pep_side(j)*side_calf(j)
+ enddo
+ cosalfa=(pep_side0pept_group/ &
+ (dist_pep_side*dist_side_calf))
+ fac_alfa_sin=1.0d0-cosalfa**2
+ fac_alfa_sin=dsqrt(fac_alfa_sin)
+ rkprim=fac_alfa_sin*(long-short)+short
+!C rkprim=short
+
+!C now costhet_grad
+ cosphi=1.0d0/dsqrt(1.0d0+rkprim**2/dist_pep_side**2)
+!C cosphi=0.6
+ cosphi_fac=cosphi**3*rkprim**2*(-0.5d0)/dist_pep_side**4
+ sinphi=rkprim/dist_pep_side/dsqrt(1.0d0+rkprim**2/ &
+ dist_pep_side**2)
+!C sinphi=0.8
+ do j=1,3
+ cosphi_grad_long(j)=cosphi_fac*pep_side(j) &
+ +cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) &
+ *(long-short)/fac_alfa_sin*cosalfa/ &
+ ((dist_pep_side*dist_side_calf))* &
+ ((side_calf(j))-cosalfa* &
+ ((pep_side(j)/dist_pep_side)*dist_side_calf))
+!C cosphi_grad_long(j)=0.0d0
+ cosphi_grad_loc(j)=cosphi**3*0.5d0/dist_pep_side**2*(-rkprim) &
+ *(long-short)/fac_alfa_sin*cosalfa &
+ /((dist_pep_side*dist_side_calf))* &
+ (pep_side(j)- &
+ cosalfa*side_calf(j)/dist_side_calf*dist_pep_side)
+!C cosphi_grad_loc(j)=0.0d0
+ enddo
+!C print *,sinphi,sinthet
+ VofOverlap=VSolvSphere/2.0d0*(1.0d0-dsqrt(1.0d0-sinphi*sinthet)) &
+ /VSolvSphere_div
+!C & *wshield
+!C now the gradient...
+ do j=1,3
+ grad_shield(j,i)=grad_shield(j,i) &
+!C gradient po skalowaniu
+ +(sh_frac_dist_grad(j)*VofOverlap &
+!C gradient po costhet
+ +scale_fac_dist*VSolvSphere/VSolvSphere_div/4.0d0* &
+ (1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*( &
+ sinphi/sinthet*costhet*costhet_grad(j) &
+ +sinthet/sinphi*cosphi*cosphi_grad_long(j))) &
+ )*wshield
+!C grad_shield_side is Cbeta sidechain gradient
+ grad_shield_side(j,ishield_list(i),i)=&
+ (sh_frac_dist_grad(j)*-2.0d0&
+ *VofOverlap&
+ -scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*&
+ (1.0d0/(-dsqrt(1.0d0-sinphi*sinthet))*(&
+ sinphi/sinthet*costhet*costhet_grad(j)&
+ +sinthet/sinphi*cosphi*cosphi_grad_long(j))) &
+ )*wshield
+! print *, 1.0d0/(-dsqrt(1.0d0-sinphi*sinthet)),&
+! sinphi/sinthet,&
+! +sinthet/sinphi,"HERE"
+ grad_shield_loc(j,ishield_list(i),i)= &
+ scale_fac_dist*VSolvSphere/VSolvSphere_div/2.0d0*&
+ (1.0d0/(dsqrt(1.0d0-sinphi*sinthet))*(&
+ sinthet/sinphi*cosphi*cosphi_grad_loc(j)&
+ ))&
+ *wshield
+! print *,grad_shield_loc(j,ishield_list(i),i)
+ enddo
+ VolumeTotal=VolumeTotal+VofOverlap*scale_fac_dist
+ enddo
+ fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield)
+
+! write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i)
+ enddo
+ return
+ end subroutine set_shield_fac2
+!----------------------------------------------------------------------------
+! SOUBROUTINE FOR AFM
+ subroutine AFMvel(Eafmforce)
+ use MD_data, only:totTafm
+ real(kind=8),dimension(3) :: diffafm,cbeg,cend
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i,j
+!C Only for check grad COMMENT if not used for checkgrad
+!C totT=3.0d0
+!C--------------------------------------------------------
+!C print *,"wchodze"
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ cbeg=0.0d0
+ cend=0.0d0
+ if (afmbeg.eq.-1) then
+ do i=1,nbegafmmat
+ do j=1,3
+ cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat
+ enddo
+ enddo
+ else
+ do j=1,3
+ cbeg(j)=c(j,afmend)
+ enddo
+ endif
+ if (afmend.eq.-1) then
+ do i=1,nendafmmat
+ do j=1,3
+ cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat
+ enddo
+ enddo
+ else
+ cend(j)=c(j,afmend)
+ endif
+
+ do i=1,3
+ diffafm(i)=cend(i)-cbeg(i)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! totTafm=3.0
+ Eafmforce=0.5d0*forceAFMconst &
+ *(distafminit+totTafm*velAFMconst-afmdist)**2
+!C Eafmforce=-forceAFMconst*(dist-distafminit)
+ if (afmend.eq.-1) then
+ do i=1,nendafmmat
+ do j=1,3
+ gradafm(j,afmendcentr(i)-1)=-forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(j)/afmdist/nendafmmat
+ enddo
+ enddo
+ else
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ enddo
+ endif
+ if (afmbeg.eq.-1) then
+ do i=1,nbegafmmat
+ do j=1,3
+ gradafm(i,afmbegcentr(i)-1)=forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ enddo
+ enddo
+ else
+ do i=1,3
+ gradafm(i,afmbeg-1)=forceAFMconst* &
+ (distafminit+totTafm*velAFMconst-afmdist) &
+ *diffafm(i)/afmdist
+ enddo
+ endif
+! print *,'AFM',Eafmforce,totTafm*velAFMconst,afmdist
+ return
+ end subroutine AFMvel
+!---------------------------------------------------------
+ subroutine AFMforce(Eafmforce)
+
+ real(kind=8),dimension(3) :: diffafm
+! real(kind=8) ::afmdist
+ real(kind=8) :: afmdist,Eafmforce
+ integer :: i
+ afmdist=0.0d0
+ Eafmforce=0.0d0
+ do i=1,3
+ diffafm(i)=c(i,afmend)-c(i,afmbeg)
+ afmdist=afmdist+diffafm(i)**2
+ enddo
+ afmdist=dsqrt(afmdist)
+! print *,afmdist,distafminit
+ Eafmforce=-forceAFMconst*(afmdist-distafminit)
+ do i=1,3
+ gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/afmdist
+ gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/afmdist
+ enddo
+!C print *,'AFM',Eafmforce
+ return
+ end subroutine AFMforce
+
+!-----------------------------------------------------------------------------
+#ifdef WHAM
+ subroutine read_ssHist
+! implicit none
+! Includes
+! include 'DIMENSIONS'
+! include "DIMENSIONS.FREE"
+! include 'COMMON.FREE'
+! Local variables
+ integer :: i,j
+ character(len=80) :: controlcard
+
+ do i=1,dyn_nssHist
+ call card_concat(controlcard,.true.)
+ read(controlcard,*) &
+ dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
+ enddo
+
+ return
+ end subroutine read_ssHist
+#endif
+!-----------------------------------------------------------------------------
+ integer function indmat(i,j)
+!el
+! get the position of the jth ijth fragment of the chain coordinate system
+! in the fromto array.
+ integer :: i,j
+
+ indmat=((2*(nres-2)-i)*(i-1))/2+j-1
+ return
+ end function indmat
+!-----------------------------------------------------------------------------
+ real(kind=8) function sigm(x)
+!el
+ real(kind=8) :: x
+ sigm=0.25d0*x
+ return
+ end function sigm
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+ subroutine alloc_ener_arrays
+!EL Allocation of arrays used by module energy
+ use MD_data, only: mset
+!el local variables
+ integer :: i,j
+
+ if(nres.lt.100) then
+ maxconts=10*nres
+ elseif(nres.lt.200) then
+ maxconts=10*nres ! Max. number of contacts per residue
+ else
+ maxconts=10*nres ! (maxconts=maxres/4)
+ endif
+ maxcont=100*nres ! Max. number of SC contacts
+ maxvar=6*nres ! Max. number of variables
+!el maxdim=(nres-1)*(nres-2)/2 ! Max. number of derivatives of virtual-bond
+ maxdim=nres*(nres-2)/2 ! Max. number of derivatives of virtual-bond
+!----------------------
+! arrays in subroutine init_int_table
+!el#ifdef MPI
+!el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1)
+!el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1)
+!el#endif
+ allocate(nint_gr(nres))
+ allocate(nscp_gr(nres))
+ allocate(ielstart(nres))
+ allocate(ielend(nres))
+!(maxres)
+ allocate(istart(nres,maxint_gr))
+ allocate(iend(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(iscpstart(nres,maxint_gr))
+ allocate(iscpend(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(ielstart_vdw(nres))
+ allocate(ielend_vdw(nres))
+!(maxres)
+ allocate(nint_gr_nucl(nres))
+ allocate(nscp_gr_nucl(nres))
+ allocate(ielstart_nucl(nres))
+ allocate(ielend_nucl(nres))
+!(maxres)
+ allocate(istart_nucl(nres,maxint_gr))
+ allocate(iend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(iscpstart_nucl(nres,maxint_gr))
+ allocate(iscpend_nucl(nres,maxint_gr))
+!(maxres,maxint_gr)
+ allocate(ielstart_vdw_nucl(nres))
+ allocate(ielend_vdw_nucl(nres))
+
+ allocate(lentyp(0:nfgtasks-1))
+!(0:maxprocs-1)
+!----------------------
+! commom.contacts
+! common /contacts/
+ if(.not.allocated(icont_ref)) allocate(icont_ref(2,maxcont))
+ allocate(icont(2,maxcont))
+!(2,maxcont)
+! common /contacts1/
+ allocate(num_cont(0:nres+4))
+!(maxres)
+#ifndef NEWCORR
+ allocate(jcont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(facont(maxconts,nres))
+!(maxconts,maxres)
+ allocate(gacont(3,maxconts,nres))
+!(3,maxconts,maxres)
+! common /contacts_hb/
+ allocate(gacontp_hb1(3,maxconts,nres))
+ allocate(gacontp_hb2(3,maxconts,nres))
+ allocate(gacontp_hb3(3,maxconts,nres))
+ allocate(gacontm_hb1(3,maxconts,nres))
+ allocate(gacontm_hb2(3,maxconts,nres))
+ allocate(gacontm_hb3(3,maxconts,nres))
+ allocate(gacont_hbr(3,maxconts,nres))
+ allocate(grij_hb_cont(3,maxconts,nres))
+ !(3,maxconts,maxres)
+ allocate(facont_hb(maxconts,nres))
+
+ allocate(ees0p(maxconts,nres))
+ allocate(ees0m(maxconts,nres))
+ allocate(d_cont(maxconts,nres))
+ allocate(ees0plist(maxconts,nres))
+
+!(maxconts,maxres)
+!(maxres)
+ allocate(jcont_hb(maxconts,nres))
+#endif
+ allocate(num_cont_hb(nres))
+!(maxconts,maxres)
+! common /rotat/
+ allocate(Ug(2,2,nres))
+ allocate(Ugder(2,2,nres))
+ allocate(Ug2(2,2,nres))
+ allocate(Ug2der(2,2,nres))
+!(2,2,maxres)
+ allocate(obrot(2,nres))
+ allocate(obrot2(2,nres))
+ allocate(obrot_der(2,nres))
+ allocate(obrot2_der(2,nres))
+!(2,maxres)
+! common /precomp1/
+ allocate(mu(2,nres))
+ allocate(muder(2,nres))
+ allocate(Ub2(2,nres))
+ Ub2(1,:)=0.0d0
+ Ub2(2,:)=0.0d0
+ allocate(Ub2der(2,nres))
+ allocate(Ctobr(2,nres))
+ allocate(Ctobrder(2,nres))
+ allocate(Dtobr2(2,nres))
+ allocate(Dtobr2der(2,nres))
+!(2,maxres)
+ allocate(EUg(2,2,nres))
+ allocate(EUgder(2,2,nres))
+ allocate(CUg(2,2,nres))
+ allocate(CUgder(2,2,nres))
+ allocate(DUg(2,2,nres))
+ allocate(Dugder(2,2,nres))
+ allocate(DtUg2(2,2,nres))
+ allocate(DtUg2der(2,2,nres))
+!(2,2,maxres)
+! common /precomp2/
+ allocate(Ug2Db1t(2,nres))
+ allocate(Ug2Db1tder(2,nres))
+ allocate(CUgb2(2,nres))
+ allocate(CUgb2der(2,nres))
+!(2,maxres)
+ allocate(EUgC(2,2,nres))
+ allocate(EUgCder(2,2,nres))
+ allocate(EUgD(2,2,nres))
+ allocate(EUgDder(2,2,nres))
+ allocate(DtUg2EUg(2,2,nres))
+ allocate(Ug2DtEUg(2,2,nres))
+!(2,2,maxres)
+ allocate(Ug2DtEUgder(2,2,2,nres))
+ allocate(DtUg2EUgder(2,2,2,nres))
+!(2,2,2,maxres)
+ allocate(b1(2,nres)) !(2,-maxtor:maxtor)
+ allocate(b2(2,nres)) !(2,-maxtor:maxtor)
+ allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
+ allocate(b2tilde(2,nres)) !(2,-maxtor:maxtor)
+
+ allocate(ctilde(2,2,nres))
+ allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
+ allocate(gtb1(2,nres))
+ allocate(gtb2(2,nres))
+ allocate(cc(2,2,nres))
+ allocate(dd(2,2,nres))
+ allocate(ee(2,2,nres))
+ allocate(gtcc(2,2,nres))
+ allocate(gtdd(2,2,nres))
+ allocate(gtee(2,2,nres))
+ allocate(gUb2(2,nres))
+ allocate(gteUg(2,2,nres))
+
+! common /rotat_old/
+ allocate(costab(nres))
+ allocate(sintab(nres))
+ allocate(costab2(nres))
+ allocate(sintab2(nres))
+!(maxres)
+! common /dipmat/
+! allocate(a_chuj(2,2,maxconts,nres))
+!(2,2,maxconts,maxres)(maxconts=maxres/4)
+! allocate(a_chuj_der(2,2,3,5,maxconts,nres))
+!(2,2,3,5,maxconts,maxres)(maxconts=maxres/4)
+! common /contdistrib/
+ allocate(ncont_sent(nres))
+ allocate(ncont_recv(nres))
+
+ allocate(iat_sent(nres))
+!(maxres)
+#ifndef NEWCORR
+ print *,"before iint_sent allocate"
+ allocate(iint_sent(4,nres,nres))
+ allocate(iint_sent_local(4,nres,nres))
+ print *,"after iint_sent allocate"
+#endif
+!(4,maxres,maxres)
+ allocate(iturn3_sent(4,0:nres+4))
+ allocate(iturn4_sent(4,0:nres+4))
+ allocate(iturn3_sent_local(4,nres))
+ allocate(iturn4_sent_local(4,nres))
+!(4,maxres)
+ allocate(itask_cont_from(0:nfgtasks-1))
+ allocate(itask_cont_to(0:nfgtasks-1))
+!(0:max_fg_procs-1)
+
+
+
+!----------------------
+! commom.deriv;
+! common /derivat/
+#ifdef NEWCORR
+ print *,"before dcdv allocate"
+ allocate(dcdv(6,nres+2))
+ allocate(dxdv(6,nres+2))
+#else
+ print *,"before dcdv allocate"
+ allocate(dcdv(6,maxdim))
+ allocate(dxdv(6,maxdim))
+#endif
+!(6,maxdim)
+ allocate(dxds(6,nres))
+!(6,maxres)
+ allocate(gradx(3,-1:nres,0:2))
+ allocate(gradc(3,-1:nres,0:2))
+!(3,maxres,2)
+ allocate(gvdwx(3,-1:nres))
+ allocate(gvdwc(3,-1:nres))
+ allocate(gelc(3,-1:nres))
+ allocate(gelc_long(3,-1:nres))
+ allocate(gvdwpp(3,-1:nres))
+ allocate(gvdwc_scpp(3,-1:nres))
+ allocate(gradx_scp(3,-1:nres))
+ allocate(gvdwc_scp(3,-1:nres))
+ allocate(ghpbx(3,-1:nres))
+ allocate(ghpbc(3,-1:nres))
+ allocate(gradcorr(3,-1:nres))
+ allocate(gradcorr_long(3,-1:nres))
+ allocate(gradcorr5_long(3,-1:nres))
+ allocate(gradcorr6_long(3,-1:nres))
+ allocate(gcorr6_turn_long(3,-1:nres))
+ allocate(gradxorr(3,-1:nres))
+ allocate(gradcorr5(3,-1:nres))
+ allocate(gradcorr6(3,-1:nres))
+ allocate(gliptran(3,-1:nres))
+ allocate(gliptranc(3,-1:nres))
+ allocate(gliptranx(3,-1:nres))
+ allocate(gshieldx(3,-1:nres))
+ allocate(gshieldc(3,-1:nres))
+ allocate(gshieldc_loc(3,-1:nres))
+ allocate(gshieldx_ec(3,-1:nres))
+ allocate(gshieldc_ec(3,-1:nres))
+ allocate(gshieldc_loc_ec(3,-1:nres))
+ allocate(gshieldx_t3(3,-1:nres))
+ allocate(gshieldc_t3(3,-1:nres))
+ allocate(gshieldc_loc_t3(3,-1:nres))
+ allocate(gshieldx_t4(3,-1:nres))
+ allocate(gshieldc_t4(3,-1:nres))
+ allocate(gshieldc_loc_t4(3,-1:nres))
+ allocate(gshieldx_ll(3,-1:nres))
+ allocate(gshieldc_ll(3,-1:nres))
+ allocate(gshieldc_loc_ll(3,-1:nres))
+ allocate(grad_shield(3,-1:nres))
+ allocate(gg_tube_sc(3,-1:nres))
+ allocate(gg_tube(3,-1:nres))
+ allocate(gradafm(3,-1:nres))
+ allocate(gradb_nucl(3,-1:nres))
+ allocate(gradbx_nucl(3,-1:nres))
+ allocate(gvdwpsb1(3,-1:nres))
+ allocate(gelpp(3,-1:nres))
+ allocate(gvdwpsb(3,-1:nres))
+ allocate(gelsbc(3,-1:nres))
+ allocate(gelsbx(3,-1:nres))
+ allocate(gvdwsbx(3,-1:nres))
+ allocate(gvdwsbc(3,-1:nres))
+ allocate(gsbloc(3,-1:nres))
+ allocate(gsblocx(3,-1:nres))
+ allocate(gradcorr_nucl(3,-1:nres))
+ allocate(gradxorr_nucl(3,-1:nres))
+ allocate(gradcorr3_nucl(3,-1:nres))
+ allocate(gradxorr3_nucl(3,-1:nres))
+ allocate(gvdwpp_nucl(3,-1:nres))
+ allocate(gradpepcat(3,-1:nres))
+ allocate(gradpepcatx(3,-1:nres))
+ allocate(gradpepmart(3,-1:nres))
+ allocate(gradpepmartx(3,-1:nres))
+ allocate(gradcatcat(3,-1:nres))
+ allocate(gradnuclcat(3,-1:nres))
+ allocate(gradnuclcatx(3,-1:nres))
+ allocate(gradlipbond(3,-1:nres))
+ allocate(gradlipang(3,-1:nres))
+ allocate(gradliplj(3,-1:nres))
+ allocate(gradlipelec(3,-1:nres))
+ allocate(gradcattranc(3,-1:nres))
+ allocate(gradcattranx(3,-1:nres))
+ allocate(gradcatangx(3,-1:nres))
+ allocate(gradcatangc(3,-1:nres))
+!(3,maxres)
+ allocate(grad_shield_side(3,maxcontsshi,-1:nres))
+ allocate(grad_shield_loc(3,maxcontsshi,-1:nres))
+! grad for shielding surroing
+ allocate(gloc(0:maxvar,0:2))
+ allocate(gloc_x(0:maxvar,2))
+!(maxvar,2)
+ allocate(gel_loc(3,-1:nres))
+ allocate(gel_loc_long(3,-1:nres))
+ allocate(gcorr3_turn(3,-1:nres))
+ allocate(gcorr4_turn(3,-1:nres))
+ allocate(gcorr6_turn(3,-1:nres))
+ allocate(gradb(3,-1:nres))
+ allocate(gradbx(3,-1:nres))
+!(3,maxres)
+ allocate(gel_loc_loc(maxvar))
+ allocate(gel_loc_turn3(maxvar))
+ allocate(gel_loc_turn4(maxvar))
+ allocate(gel_loc_turn6(maxvar))
+ allocate(gcorr_loc(maxvar))
+ allocate(g_corr5_loc(maxvar))
+ allocate(g_corr6_loc(maxvar))
+!(maxvar)
+ allocate(gsccorc(3,-1:nres))
+ allocate(gsccorx(3,-1:nres))
+!(3,maxres)
+ allocate(gsccor_loc(-1:nres))
+!(maxres)
+ allocate(gvdwx_scbase(3,-1:nres))
+ allocate(gvdwc_scbase(3,-1:nres))
+ allocate(gvdwx_pepbase(3,-1:nres))
+ allocate(gvdwc_pepbase(3,-1:nres))
+ allocate(gvdwx_scpho(3,-1:nres))
+ allocate(gvdwc_scpho(3,-1:nres))
+ allocate(gvdwc_peppho(3,-1:nres))
+
+ allocate(dtheta(3,2,-1:nres))
+!(3,2,maxres)
+ allocate(gscloc(3,-1:nres))
+ allocate(gsclocx(3,-1:nres))
+!(3,maxres)
+ allocate(dphi(3,3,-1:nres))
+ allocate(dalpha(3,3,-1:nres))
+ allocate(domega(3,3,-1:nres))
+!(3,3,maxres)
+! common /deriv_scloc/
+ allocate(dXX_C1tab(3,nres))
+ allocate(dYY_C1tab(3,nres))
+ allocate(dZZ_C1tab(3,nres))
+ allocate(dXX_Ctab(3,nres))
+ allocate(dYY_Ctab(3,nres))
+ allocate(dZZ_Ctab(3,nres))
+ allocate(dXX_XYZtab(3,nres))
+ allocate(dYY_XYZtab(3,nres))
+ allocate(dZZ_XYZtab(3,nres))
+!(3,maxres)
+! common /mpgrad/
+ allocate(jgrad_start(nres))
+ allocate(jgrad_end(nres))
+!(maxres)
+!----------------------
+
+! common /indices/
+ allocate(ibond_displ(0:nfgtasks-1))
+ allocate(ibond_count(0:nfgtasks-1))
+ allocate(ithet_displ(0:nfgtasks-1))
+ allocate(ithet_count(0:nfgtasks-1))
+ allocate(iphi_displ(0:nfgtasks-1))
+ allocate(iphi_count(0:nfgtasks-1))
+ allocate(iphi1_displ(0:nfgtasks-1))
+ allocate(iphi1_count(0:nfgtasks-1))
+ allocate(ivec_displ(0:nfgtasks-1))
+ allocate(ivec_count(0:nfgtasks-1))
+ allocate(iset_displ(0:nfgtasks-1))
+ allocate(iset_count(0:nfgtasks-1))
+ allocate(iint_count(0:nfgtasks-1))
+ allocate(iint_displ(0:nfgtasks-1))
+!(0:max_fg_procs-1)
+!----------------------
+! common.MD
+! common /mdgrad/
+ allocate(gcart(3,-1:nres))
+ allocate(gxcart(3,-1:nres))
+!(3,0:MAXRES)
+ allocate(gradcag(3,-1:nres))
+ allocate(gradxag(3,-1:nres))
+!(3,MAXRES)
+! common /back_constr/
+!el in energy:Econstr_back allocate((:),allocatable :: utheta,ugamma,uscdiff !(maxfrag_back)
+ allocate(dutheta(nres))
+ allocate(dugamma(nres))
+!(maxres)
+ allocate(duscdiff(3,-1:nres))
+ allocate(duscdiffx(3,-1:nres))
+!(3,maxres)
+!el i io:read_fragments
+! allocate((:,:,:),allocatable :: wfrag_back !(3,maxfrag_back,maxprocs/20)
+! allocate((:,:,:),allocatable :: ifrag_back !(3,maxfrag_back,maxprocs/20)
+! common /qmeas/
+! allocate(qinfrag(50,nprocs/20),wfrag(50,nprocs/20)) !(50,maxprocs/20)
+! allocate(qinpair(100,nprocs/20),wpair(100,nprocs/20)) !(100,maxprocs/20)
+ allocate(mset(0:nprocs)) !(maxprocs/20)
+ mset(:)=0
+! allocate(ifrag(2,50,nprocs/20)) !(2,50,maxprocs/20)
+! allocate(ipair(2,100,nprocs/20)) !(2,100,maxprocs/20)
+ allocate(dUdconst(3,0:nres))
+ allocate(dUdxconst(3,0:nres))
+ allocate(dqwol(3,0:nres))
+ allocate(dxqwol(3,0:nres))
+!(3,0:MAXRES)
+!----------------------
+! common.sbridge
+! common /sbridge/ in io_common: read_bridge
+!el allocate((:),allocatable :: iss !(maxss)
+! common /links/ in io_common: read_bridge
+!el real(kind=8),dimension(:),allocatable :: dhpb,forcon,dhpb1 !(maxdim) !el dhpb1 !!! nie używane
+!el integer,dimension(:),allocatable :: ihpb,jhpb,ibecarb !(maxdim) !el ibecarb !!! nie używane
+! common /dyn_ssbond/
+! and side-chain vectors in theta or phi.
+ allocate(dyn_ssbond_ij(10000))
+!(maxres,maxres)
+! do i=1,nres
+! do j=i+1,nres
+ dyn_ssbond_ij(:)=1.0d300
+! enddo
+! enddo
+
+! if (nss.gt.0) then
+ allocate(idssb(maxdim),jdssb(maxdim))
+! allocate(newihpb(nss),newjhpb(nss))
+!(maxdim)
+! endif
+ allocate(ishield_list(-1:nres))
+ allocate(shield_list(maxcontsshi,-1:nres))
+ allocate(dyn_ss_mask(nres))
+ allocate(fac_shield(-1:nres))
+ allocate(enetube(nres*2))
+ allocate(enecavtube(nres*2))
+
+!(maxres)
+ dyn_ss_mask(:)=.false.
+!----------------------
+! common.sccor
+! Parameters of the SCCOR term
+! common/sccor/
+!el in io_conf: parmread
+! allocate(v1sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp))
+! allocate(v2sccor(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
+! allocate(v0sccor(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
+! allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
+! allocate(nterm_sccor(-ntyp:ntyp,-ntyp:ntyp))
+! allocate(nlor_sccor(-ntyp:ntyp,-ntyp:ntyp)) !(-ntyp:ntyp,-ntyp:ntyp)
+! allocate(vlor1sccor(maxterm_sccor,20,20))
+! allocate(vlor2sccor(maxterm_sccor,20,20))
+! allocate(vlor3sccor(maxterm_sccor,20,20)) !(maxterm_sccor,20,20)
+!----------------
+ allocate(gloc_sc(3,0:2*nres,0:10))
+!(3,0:maxres2,10)maxres2=2*maxres
+ allocate(dcostau(3,3,3,2*nres))
+ allocate(dsintau(3,3,3,2*nres))
+ allocate(dtauangle(3,3,3,2*nres))
+ allocate(dcosomicron(3,3,3,2*nres))
+ allocate(domicron(3,3,3,2*nres))
+!(3,3,3,maxres2)maxres2=2*maxres
+!----------------------
+! common.var
+! common /restr/
+ allocate(varall(maxvar))
+!(maxvar)(maxvar=6*maxres)
+ allocate(mask_theta(nres))
+ allocate(mask_phi(nres))
+ allocate(mask_side(nres))
+!(maxres)
+!----------------------
+! common.vectors
+! common /vectors/
+ allocate(uy(3,nres))
+ allocate(uz(3,nres))
+!(3,maxres)
+ allocate(uygrad(3,3,2,nres))
+ allocate(uzgrad(3,3,2,nres))
+!(3,3,2,maxres)
+ print *,"before all 300"
+! allocateion of lists JPRDLA
+ allocate(newcontlistppi(300*nres))
+ allocate(newcontlistscpi(350*nres))
+ allocate(newcontlisti(300*nres))
+ allocate(newcontlistppj(300*nres))
+ allocate(newcontlistscpj(350*nres))
+ allocate(newcontlistj(300*nres))
+ allocate(newcontlistmartpi(300*nres))
+ allocate(newcontlistmartpj(300*nres))
+ allocate(newcontlistmartsci(300*nres))
+ allocate(newcontlistmartscj(300*nres))
+
+ allocate(newcontlistcatsctrani(300*nres))
+ allocate(newcontlistcatsctranj(300*nres))
+ allocate(newcontlistcatptrani(300*nres))
+ allocate(newcontlistcatptranj(300*nres))
+ allocate(newcontlistcatscnormi(300*nres))
+ allocate(newcontlistcatscnormj(300*nres))
+ allocate(newcontlistcatpnormi(300*nres))
+ allocate(newcontlistcatpnormj(300*nres))
+ allocate(newcontlistcatcatnormi(900*nres))
+ allocate(newcontlistcatcatnormj(900*nres))
+
+ allocate(newcontlistcatscangi(300*nres))
+ allocate(newcontlistcatscangj(300*nres))
+ allocate(newcontlistcatscangfi(300*nres))
+ allocate(newcontlistcatscangfj(300*nres))
+ allocate(newcontlistcatscangfk(300*nres))
+ allocate(newcontlistcatscangti(300*nres))
+ allocate(newcontlistcatscangtj(300*nres))
+ allocate(newcontlistcatscangtk(300*nres))
+ allocate(newcontlistcatscangtl(300*nres))
+
+
+ return
+ end subroutine alloc_ener_arrays
+!-----------------------------------------------------------------
+ subroutine ebond_nucl(estr_nucl)
+!c
+!c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds
+!c
+
+ real(kind=8),dimension(3) :: u,ud
+ real(kind=8) :: usum,uprod,uprod1,uprod2,usumsqder
+ real(kind=8) :: estr_nucl,diff
+ integer :: iti,i,j,k,nbi
+ estr_nucl=0.0d0
+!C print *,"I enter ebond"
+ if (energy_dec) &
+ write (iout,*) "ibondp_start,ibondp_end",&
+ ibondp_nucl_start,ibondp_nucl_end
+ do i=ibondp_nucl_start,ibondp_nucl_end
+
+ if (itype(i-1,2).eq.ntyp1_molec(2)&
+ .and.itype(i,2).eq.ntyp1_molec(2)) cycle
+ if (itype(i-1,2).eq.ntyp1_molec(2)&
+ .or. itype(i,2).eq.ntyp1_molec(2)) then
+!C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+!C do j=1,3
+!C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
+!C *dc(j,i-1)/vbld(i)
+!C enddo
+!C if (energy_dec) write(iout,*) &
+!C "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+ diff = vbld(i)-vbldpDUM
+ else
+ diff = vbld(i)-vbldp0_nucl
projects / unres4.git / blobdiff