PROGRAM chpsi C empirical correction to vacuum peptide FF c on force field belonging to atoms 1...N4 c based on potential on atoms 1..N6 c nat0 - maximum number of atoms c m0 - max number of clusters c nd0 - max dimension of the inversion problem IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (nat0=800,MX3=3*nat0,MX=nat0,nd0=10) PARAMETER (nat0big=13000) DIMENSION r(3*nat0big),ity(nat0big),ind(nat0big), 2rv(MX3),ityv(nat0),i4(nat0),fi(nat0) logical lex,lper,ltinker character*80 ts,name1 character*2 at DATA IERR,n6/0,6/ DATA qh,qo/0.42d0,-0.84d0/ C WRITE(6,60000) 60000 FORMAT(' SOLVENT ELSTAT FIELD ',/,/, 1 ' WITH SIDECHAIN CHARGES ',/, 1 ' (NH3+, COO- ONLY) ',/,/, 1 ' Input: SMALL.X ',/, 1 ' CCT.INP (atom map) ',/, 1 ' HYPEP.OPT (optional) ',/, 1 ' Output: FILE.PSI ',/,/) C OPEN(16,file='FILE.PSI') c OPEN(20,FILE='SMALL.X',STATUS='OLD') read(20,*) read(20,*)natv if(natv.gt.nat0)call report('Too many atoms') do 41 ia=1,natv 41 read(20,*)ityv(ia),(rv(3*(ia-1)+ix),ix=1,3) close(20) c lper=.false. ltinker=.false. name1='BIG.X' inquire(file='HYPEP.OPT',exist=lex) if(lex)then open(55,file='HYPEP.OPT') 1001 read(55,2000,end=1002,err=1002)ts 2000 format(a80) if(ts(1:4).eq.'NAME')read(55,2000)name1 if(ts(1:4).eq.'TINK')then ltinker=.true. write(6,*)'Tinker format for BIG' name1='tinker.xyz' endif if(ts(1:4).eq.'PERI')then lper=.true. read(55,*)DPCX,DPCY,DPCZ dx2=DPCX/2.0d0 dy2=DPCY/2.0d0 dz2=DPCZ/2.0d0 write(6,*)' Periodic boundary conditions' endif goto 1001 1002 close(55) write(6,*)' HYPEP.OPT found' endif c OPEN(20,FILE=name1,STATUS='OLD') if(.not.ltinker)read(20,*) read(20,*)nat if(nat.gt.nat0big)call report('Too many atoms') if(ltinker)then do 4141 ia=1,nat read(20,2020)at,(r(3*(ia-1)+ix),ix=1,3) 2020 format(6x,2x,a2,1x,3f12.6) ity(ia)=0 if(at(1:1).eq.'H')ity(ia)=1 if(at(1:1).eq.'C')ity(ia)=6 if(at(1:1).eq.'N')then ity(ia)=7 if(at(2:2).eq.'2')then ity(ia)=72 else if(at(2:2).eq.'3')then ity(ia)=73 endif endif endif if(at(1:1).eq.'O')then ity(ia)=8 if(at(2:2).eq.'2')then ity(ia)=82 endif endif 4141 if(ity(ia).eq.0)write(6,*)'Unknown atom ',ia,at else do 414 ia=1,nat 414 read(20,*)ity(ia),(r(3*(ia-1)+ix),ix=1,3) endif close(20) write(6,*)'Geometry read for BIG.' c OPEN(20,FILE='CCT.INP',STATUS='OLD') read(20,*) read(20,*) read(20,201)(ind(i),i=1,nat) read(20,201)(ind(i),i=1,nat) 201 format(20i3) close(20) write(6,*)' Atom map read from CCT.INP' c c identify amide groups in peptide iag=0 do 1 ia=1,nat if(ity(ia).eq.6)then xi=r(3*(ia-1)+1) yi=r(3*(ia-1)+2) zi=r(3*(ia-1)+3) iin=0 iih=0 iian=0 iiac=0 iio=0 in=0 ih=0 ian=0 iac=0 io=0 it=0 do 2 ja=1,nat if(ja.ne.ia)then xj=r(3*(ja-1)+1) yj=r(3*(ja-1)+2) zj=r(3*(ja-1)+3) d2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 if(d2.lt.2.7d0)then it=it+1 if(ity(ja).eq.8)then io=io+1 iio=ja endif if(ity(ja).eq.6)then iac=iac+1 iiac=ja endif if(ity(ja).eq.7)then in=in+1 iin=ja endif endif endif 2 continue if(it.eq.3.and.io.eq.1.and.iac.eq.1.and.in.eq.1)then xn=r(3*(iin-1)+1) yn=r(3*(iin-1)+2) zn=r(3*(iin-1)+3) do 3 ja=1,nat if(ja.ne.iin)then xj=r(3*(ja-1)+1) yj=r(3*(ja-1)+2) zj=r(3*(ja-1)+3) d2=(xn-xj)**2+(yn-yj)**2+(zn-zj)**2 if(d2.lt.2.7d0)then if(ity(ja).eq.1)then ih=ih+1 iih=ja endif if(ity(ja).eq.6.and.ja.ne.ia)then ian=ian+1 iian=ja endif endif endif 3 continue if(ih.eq.1.and.ian.eq.1)then iag=iag+1 iic=ia write(6,600)iag,iih,iin,iic,iio,iian,iiac write(6,601)ind(iih),ind(iin),ind(iic),ind(iio),ind(iian), 1 ind(iiac) 600 format(' Amide group number ',i3,': ',6i3) 601 format(' Vacuum system ',3x,': ',6i3) if(ind(iih).eq.0.or.ind(iin).eq.0.or.ind(iic).eq.0.or. 1 ind(iio).eq.0.or.ind(iian).eq.0.or.ind(iiac).eq.0)then write(6,65) 65 format('Amide group ill-mapped') else i4(1)=iih i4(2)=iin i4(3)=iic i4(4)=iio if(n4.eq.6)then i4(5)=iian i4(6)=iiac endif c c find waters and their potentials do 4 i=1,n6 4 fi(i)=0.0d0 do 6 jj=1,nat if(ity(jj).eq.8)then j1=3*(jj-1) xo=r(j1+1) yo=r(j1+2) zo=r(j1+3) nh=0 do 7 kk=1,nat if(ity(kk).eq.1)then k1=3*(kk-1) xh=r(k1+1) yh=r(k1+2) zh=r(k1+3) d2=(xo-xh)**2+(yo-yh)**2+(zo-zh)**2 if(d2.gt.0.64d0.and.d2.lt.1.44d0)then nh=nh+1 if(nh.gt.2)call report('Too many Hs around one O') if(nh.eq.1)then x1=xh y1=yh z1=zh else x2=xh y2=yh z2=zh endif endif endif 7 continue if(nh.eq.2)then call add(iih ,1,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iin ,2,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iic ,3,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iio ,4,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iian,5,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iiac,6,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) endif endif 6 continue do 66 jj=1,nat if(ity(jj).eq.72.or.ity(jj).eq.73.or.ity(jj).eq.82)then j1=3*(jj-1) x=r(j1+1) y=r(j1+2) z=r(j1+3) if(ity(jj).eq.72)qsc=0.5d0 if(ity(jj).eq.73)qsc=1.0d0 if(ity(jj).eq.82)qsc=-0.5d0 call add(iih ,1,r,x,y,z,x,y,z,x,y,z,fi,0.0d0,qsc, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iin ,2,r,x,y,z,x,y,z,x,y,z,fi,0.0d0,qsc, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iic ,3,r,x,y,z,x,y,z,x,y,z,fi,0.0d0,qsc, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iio ,4,r,x,y,z,x,y,z,x,y,z,fi,0.0d0,qsc, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iian,5,r,x,y,z,x,y,z,x,y,z,fi,0.0d0,qsc, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) call add(iiac,6,r,x,y,z,x,y,z,x,y,z,fi,0.0d0,qsc, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) endif 66 continue write(6,1603)(fi(ii),ii=1,N6) write(16,1605)(fi(ii),ii=1,N6) 1603 format(8f9.3) 1605 format(6f12.6) c c endif endif endif endif 1 continue write(6,*)iag,' amide groups found' c close(16) stop c end subroutine report(s) character*(*) s write(6,*)s stop end subroutine add(ii,ic,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo, 1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ) real*8 r(*),xo,yo,zo,x1,y1,z1,x2,y2,z2,fi(*),x,y,z,ro,r1,r2,qh,qo, 1dx2,dy2,dz2,DPCX,DPCY,DPCZ integer*4 ii,ic logical lper c x=r(3*(ii-1)+1) y=r(3*(ii-1)+2) z=r(3*(ii-1)+3) if(lper)then if(x-xo.gt.dx2.and.x-x1.gt.dx2.and.x-x2.gt.dx2)x=x-DPCX if(y-yo.gt.dy2.and.y-y1.gt.dy2.and.y-y2.gt.dy2)y=y-DPCY if(z-zo.gt.dz2.and.z-z1.gt.dz2.and.z-z2.gt.dz2)z=z-DPCZ if(xo-x.gt.dx2.and.x1-x.gt.dx2.and.x2-x.gt.dx2)x=x+DPCX if(yo-y.gt.dy2.and.y1-y.gt.dy2.and.y2-y.gt.dy2)y=y+DPCY if(zo-z.gt.dz2.and.z1-z.gt.dz2.and.z2-z.gt.dz2)z=z+DPCZ endif ro=dsqrt((xo-x)**2+(yo-y)**2+(zo-z)**2) r1=dsqrt((x1-x)**2+(y1-y)**2+(z1-z)**2) r2=dsqrt((x2-x)**2+(y2-y)**2+(z2-z)**2) c fi(ic)=fi(ic)+qh/r1+qh/r2+qo/ro return end