--- /dev/null
+ subroutine eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4)
+C
+C This subroutine calculates the average interaction energy and its gradient
+C in the virtual-bond vectors between non-adjacent peptide groups, based on
+C the potential described in Liwo et al., Protein Sci., 1993, 2, 1715.
+C The potential depends both on the distance of peptide-group centers and on
+C the orientation of the CA-CA virtual bonds.
+C
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.TORSION'
+ include 'COMMON.VECTORS'
+ include 'COMMON.FFIELD'
+ dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+ & erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
+ double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
+ & aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+ common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
+c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
+#ifdef MOMENT
+ double precision scal_el /1.0d0/
+#else
+ double precision scal_el /0.5d0/
+#endif
+C 12/13/98
+C 13-go grudnia roku pamietnego...
+ double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
+ & 0.0d0,1.0d0,0.0d0,
+ & 0.0d0,0.0d0,1.0d0/
+cd write(iout,*) 'In EELEC'
+cd do i=1,nloctyp
+cd write(iout,*) 'Type',i
+cd write(iout,*) 'B1',B1(:,i)
+cd write(iout,*) 'B2',B2(:,i)
+cd write(iout,*) 'CC',CC(:,:,i)
+cd write(iout,*) 'DD',DD(:,:,i)
+cd write(iout,*) 'EE',EE(:,:,i)
+cd enddo
+cd call check_vecgrad
+cd stop
+ if (icheckgrad.eq.1) then
+ do i=1,nres-1
+ fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
+ do k=1,3
+ dc_norm(k,i)=dc(k,i)*fac
+ enddo
+c write (iout,*) 'i',i,' fac',fac
+ enddo
+ endif
+ if (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
+ & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 .or.
+ & wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) then
+c call vec_and_deriv
+ call set_matrices
+ endif
+cd do i=1,nres-1
+cd write (iout,*) 'i=',i
+cd do k=1,3
+cd write (iout,'(i5,2f10.5)') k,uy(k,i),uz(k,i)
+cd enddo
+cd do k=1,3
+cd write (iout,'(f10.5,2x,3f10.5,2x,3f10.5)')
+cd & uz(k,i),(uzgrad(k,l,1,i),l=1,3),(uzgrad(k,l,2,i),l=1,3)
+cd enddo
+cd enddo
+ num_conti_hb=0
+ ees=0.0D0
+ evdw1=0.0D0
+ eel_loc=0.0d0
+ eello_turn3=0.0d0
+ eello_turn4=0.0d0
+ ind=0
+ do i=1,nres
+ num_cont_hb(i)=0
+ enddo
+cd print '(a)','Enter EELEC'
+cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
+ do i=1,nres
+ gel_loc_loc(i)=0.0d0
+ gcorr_loc(i)=0.0d0
+ enddo
+ do i=iatel_s,iatel_e
+ dxi=dc(1,i)
+ dyi=dc(2,i)
+ dzi=dc(3,i)
+ dx_normi=dc_norm(1,i)
+ dy_normi=dc_norm(2,i)
+ dz_normi=dc_norm(3,i)
+ xmedi=c(1,i)+0.5d0*dxi
+ ymedi=c(2,i)+0.5d0*dyi
+ zmedi=c(3,i)+0.5d0*dzi
+ num_conti=0
+c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+ do j=ielstart(i),ielend(i)
+ ind=ind+1
+ iteli=itel(i)
+ itelj=itel(j)
+ if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+ aaa=app(iteli,itelj)
+ bbb=bpp(iteli,itelj)
+ ael6i=ael6(iteli,itelj)
+ ael3i=ael3(iteli,itelj)
+C Diagnostics only!!!
+c aaa=0.0D0
+c bbb=0.0D0
+c ael6i=0.0D0
+c ael3i=0.0D0
+C End diagnostics
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+ xj=c(1,j)+0.5D0*dxj-xmedi
+ yj=c(2,j)+0.5D0*dyj-ymedi
+ zj=c(3,j)+0.5D0*dzj-zmedi
+ rij=xj*xj+yj*yj+zj*zj
+ rrmij=1.0D0/rij
+ rij=dsqrt(rij)
+ rmij=1.0D0/rij
+ r3ij=rrmij*rmij
+ r6ij=r3ij*r3ij
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ fac=cosa-3.0D0*cosb*cosg
+ ev1=aaa*r6ij*r6ij
+c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
+ if (j.eq.i+2) ev1=scal_el*ev1
+ ev2=bbb*r6ij
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ evdwij=ev1+ev2
+ el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
+ el2=fac4*fac
+ eesij=el1+el2
+C 12/26/95 - for the evaluation of multi-body H-bonding interactions
+ ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
+ ees=ees+eesij
+ evdw1=evdw1+evdwij
+cd write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
+cd & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
+cd & 1.0D0/dsqrt(rrmij),evdwij,eesij,
+cd & xmedi,ymedi,zmedi,xj,yj,zj
+
+ if (energy_dec) then
+ write (iout,'(a6,2i5,0pf7.3)') 'evdw1',i,j,evdwij
+ write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+ endif
+
+C
+C Calculate contributions to the Cartesian gradient.
+C
+#ifdef SPLITELE
+ facvdw=-6*rrmij*(ev1+evdwij)
+ facel=-3*rrmij*(el1+eesij)
+ fac1=fac
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
+*
+* Radial derivatives. First process both termini of the fragment (i,j)
+*
+ ggg(1)=facel*xj
+ ggg(2)=facel*yj
+ ggg(3)=facel*zj
+ do k=1,3
+ ghalf=0.5D0*ggg(k)
+ gelc(k,i)=gelc(k,i)+ghalf
+ gelc(k,j)=gelc(k,j)+ghalf
+ enddo
+*
+* Loop over residues i+1 thru j-1.
+*
+caug8 do k=i+1,j-1
+caug8 do l=1,3
+caug8 gelc(l,k)=gelc(l,k)+ggg(l)
+caug8 enddo
+caug8 enddo
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
+ do k=1,3
+ ghalf=0.5D0*ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)+ghalf
+ gvdwpp(k,j)=gvdwpp(k,j)+ghalf
+ enddo
+*
+* Loop over residues i+1 thru j-1.
+*
+caug8 do k=i+1,j-1
+caug8 do l=1,3
+caug8 gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
+caug8 enddo
+caug8 enddo
+#else
+ facvdw=ev1+evdwij
+ facel=el1+eesij
+ fac1=fac
+ fac=-3*rrmij*(facvdw+facvdw+facel)
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
+*
+* Radial derivatives. First process both termini of the fragment (i,j)
+*
+ ggg(1)=fac*xj
+ ggg(2)=fac*yj
+ ggg(3)=fac*zj
+ do k=1,3
+ ghalf=0.5D0*ggg(k)
+ gelc(k,i)=gelc(k,i)+ghalf
+ gelc(k,j)=gelc(k,j)+ghalf
+ enddo
+*
+* Loop over residues i+1 thru j-1.
+*
+caug8 do k=i+1,j-1
+caug8 do l=1,3
+caug8 gelc(l,k)=gelc(l,k)+ggg(l)
+caug8 enddo
+caug8 enddo
+#endif
+*
+* Angular part
+*
+ ecosa=2.0D0*fac3*fac1+fac4
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
+ ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
+ do k=1,3
+ dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
+ dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
+ enddo
+cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
+cd & (dcosg(k),k=1,3)
+ do k=1,3
+ ggg(k)=ecosb*dcosb(k)+ecosg*dcosg(k)
+ enddo
+ do k=1,3
+ ghalf=0.5D0*ggg(k)
+ gelc(k,i)=gelc(k,i)+ghalf
+ & +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ gelc(k,j)=gelc(k,j)+ghalf
+ & +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ enddo
+caug8 do k=i+1,j-1
+caug8 do l=1,3
+caug8 gelc(l,k)=gelc(l,k)+ggg(l)
+caug8 enddo
+caug8 enddo
+
+ enddo ! j
+ num_cont_hb(i)=num_conti
+ enddo ! i
+c write (iout,*) "Number of loop steps in EELEC:",ind
+cd do i=1,nres
+cd write (iout,'(i3,3f10.5,5x,3f10.5)')
+cd & i,(gel_loc(k,i),k=1,3),gel_loc_loc(i)
+cd enddo
+c 12/7/99 Adam eello_turn3 will be considered as a separate energy term
+ccc eel_loc=eel_loc+eello_turn3
+ return
+ end
+
projects / unres.git / blobdiff