+++ /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