--- /dev/null
+ subroutine elecont(lprint,ncont,icont,ist,ien)
+ implicit none
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.COMPAR'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.NAMES'
+ include 'COMMON.LOCAL'
+ logical lprint
+ integer i,j,k,ist,ien,iteli,itelj,ind,i1,i2,it1,it2,ic1,ic2
+ double precision rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,
+ & xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,
+ & vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,
+ & eesij,ees,evdw,ene
+ double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
+ & appc(2,2),bppc(2,2)
+ double precision elcutoff,elecutoff_14
+ integer ncont,icont(2,maxcont)
+ double precision econt(maxcont)
+*
+* Load the constants of peptide bond - peptide bond interactions.
+* Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
+* proline) - determined by averaging ECEPP energy.
+*
+* as of 7/06/91.
+*
+c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+c data rpp / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
+ data elpp6c /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
+ data elpp3c / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
+ data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
+ ees=0.0d0
+ evdw=0.0d0
+ if (lprint) write (iout,'(a)')
+ & "Constants of electrostatic interaction energy expression."
+ do i=1,2
+ do j=1,2
+ rri=rpp(i,j)**6
+ appc(i,j)=epp(i,j)*rri*rri
+ bppc(i,j)=-2.0*epp(i,j)*rri
+ ael6c(i,j)=elpp6c(i,j)*4.2**6
+ ael3c(i,j)=elpp3c(i,j)*4.2**3
+ if (lprint)
+ & write (iout,'(2i2,4e15.4)') i,j,appc(i,j),bppc(i,j),ael6c(i,j),
+ & ael3c(i,j)
+ enddo
+ enddo
+ ncont=0
+ do 1 i=ist,ien-2
+ xi=c(1,i)
+ yi=c(2,i)
+ zi=c(3,i)
+ dxi=c(1,i+1)-c(1,i)
+ dyi=c(2,i+1)-c(2,i)
+ dzi=c(3,i+1)-c(3,i)
+ xmedi=xi+0.5*dxi
+ ymedi=yi+0.5*dyi
+ zmedi=zi+0.5*dzi
+ do 4 j=i+2,ien-1
+ ind=ind+1
+ iteli=itel(i)
+ itelj=itel(j)
+ if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+ if (iteli.eq.2 .and. itelj.eq.2) goto 4
+ aaa=appc(iteli,itelj)
+ bbb=bppc(iteli,itelj)
+ ael6i=ael6c(iteli,itelj)
+ ael3i=ael3c(iteli,itelj)
+ dxj=c(1,j+1)-c(1,j)
+ dyj=c(2,j+1)-c(2,j)
+ dzj=c(3,j+1)-c(3,j)
+ xj=c(1,j)+0.5*dxj-xmedi
+ yj=c(2,j)+0.5*dyj-ymedi
+ zj=c(3,j)+0.5*dzj-zmedi
+ rrmij=1.0/(xj*xj+yj*yj+zj*zj)
+ rmij=sqrt(rrmij)
+ r3ij=rrmij*rmij
+ r6ij=r3ij*r3ij
+ vrmij=vblinv*rmij
+ cosa=(dxi*dxj+dyi*dyj+dzi*dzj)*vblinv2
+ cosb=(xj*dxi+yj*dyi+zj*dzi)*vrmij
+ cosg=(xj*dxj+yj*dyj+zj*dzj)*vrmij
+ fac=cosa-3.0*cosb*cosg
+ ev1=aaa*r6ij*r6ij
+ ev2=bbb*r6ij
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ evdwij=ev1+ev2
+ el1=fac3*(4.0+fac*fac-3.0*(cosb*cosb+cosg*cosg))
+ el2=fac4*fac
+ eesij=el1+el2
+ if (j.gt.i+2 .and. eesij.le.elcutoff .or.
+ & j.eq.i+2 .and. eesij.le.elecutoff_14) then
+ ncont=ncont+1
+ icont(1,ncont)=i
+ icont(2,ncont)=j
+ econt(ncont)=eesij
+ endif
+ ees=ees+eesij
+ evdw=evdw+evdwij
+ 4 continue
+ 1 continue
+ if (lprint) then
+ write (iout,*) 'Total average electrostatic energy: ',ees
+ write (iout,*) 'VDW energy between peptide-group centers: ',evdw
+ write (iout,*)
+ write (iout,*) 'Electrostatic contacts before pruning: '
+ do i=1,ncont
+ i1=icont(1,i)
+ i2=icont(2,i)
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+ & i,restyp(it1),i1,restyp(it2),i2,econt(i)
+ enddo
+ endif
+c For given residues keep only the contacts with the greatest energy.
+ i=0
+ do while (i.lt.ncont)
+ i=i+1
+ ene=econt(i)
+ ic1=icont(1,i)
+ ic2=icont(2,i)
+ j=i
+ do while (j.lt.ncont)
+ j=j+1
+ if (ic1.eq.icont(1,j).and.iabs(icont(2,j)-ic2).le.2 .or.
+ & ic2.eq.icont(2,j).and.iabs(icont(1,j)-ic1).le.2) then
+c write (iout,*) "i",i," j",j," ic1",ic1," ic2",ic2,
+c & " jc1",icont(1,j)," jc2",icont(2,j)," ncont",ncont
+ if (econt(j).lt.ene .and. icont(2,j).ne.icont(1,j)+2) then
+ if (ic1.eq.icont(1,j)) then
+ do k=1,ncont
+ if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.icont(2,j)
+ & .and. iabs(icont(1,k)-ic1).le.2 .and.
+ & econt(k).lt.econt(j) ) goto 21
+ enddo
+ else if (ic2.eq.icont(2,j) ) then
+ do k=1,ncont
+ if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.icont(1,j)
+ & .and. iabs(icont(2,k)-ic2).le.2 .and.
+ & econt(k).lt.econt(j) ) goto 21
+ enddo
+ endif
+c Remove ith contact
+ do k=i+1,ncont
+ icont(1,k-1)=icont(1,k)
+ icont(2,k-1)=icont(2,k)
+ econt(k-1)=econt(k)
+ enddo
+ i=i-1
+ ncont=ncont-1
+c write (iout,*) "ncont",ncont
+c do k=1,ncont
+c write (iout,*) icont(1,k),icont(2,k)
+c enddo
+ goto 20
+ else if (econt(j).gt.ene .and. ic2.ne.ic1+2)
+ & then
+ if (ic1.eq.icont(1,j)) then
+ do k=1,ncont
+ if (k.ne.i .and. k.ne.j .and. icont(2,k).eq.ic2
+ & .and. iabs(icont(1,k)-icont(1,j)).le.2 .and.
+ & econt(k).lt.econt(i) ) goto 21
+ enddo
+ else if (ic2.eq.icont(2,j) ) then
+ do k=1,ncont
+ if (k.ne.i .and. k.ne.j .and. icont(1,k).eq.ic1
+ & .and. iabs(icont(2,k)-icont(2,j)).le.2 .and.
+ & econt(k).lt.econt(i) ) goto 21
+ enddo
+ endif
+c Remove jth contact
+ do k=j+1,ncont
+ icont(1,k-1)=icont(1,k)
+ icont(2,k-1)=icont(2,k)
+ econt(k-1)=econt(k)
+ enddo
+ ncont=ncont-1
+c write (iout,*) "ncont",ncont
+c do k=1,ncont
+c write (iout,*) icont(1,k),icont(2,k)
+c enddo
+ j=j-1
+ endif
+ endif
+ 21 continue
+ enddo
+ 20 continue
+ enddo
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) 'Electrostatic contacts after pruning: '
+ do i=1,ncont
+ i1=icont(1,i)
+ i2=icont(2,i)
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(i3,2x,a,i4,2x,a,i4,f10.5)')
+ & i,restyp(it1),i1,restyp(it2),i2,econt(i)
+ enddo
+ endif
+ return
+ end
projects / unres.git / blobdiff