X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Felecont.f;fp=source%2Fwham%2Fsrc-NEWSC%2Felecont.f;h=1eff2f126c74ac2a33732d4f90be99d1e1dcbecb;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/wham/src-NEWSC/elecont.f b/source/wham/src-NEWSC/elecont.f new file mode 100755 index 0000000..1eff2f1 --- /dev/null +++ b/source/wham/src-NEWSC/elecont.f @@ -0,0 +1,207 @@ + 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