C Maximum number of terms in SC bond-stretching potential
integer maxbondterm
parameter (maxbondterm=3)
+C Maximum number of generated conformations
+ integer mxio
+ parameter (mxio=2)
#INSTALL_DIR = /opt/cray/mpt/7.3.2/gni/mpich-intel/15.0
FC = ftn
-OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic
-#OPT = -CB -g -mcmodel=medium -shared-intel -dynamic
+#OPT = -O3 -ip -mcmodel=medium -shared-intel -dynamic
+OPT = -CB -g -mcmodel=medium -shared-intel -dynamic
FFLAGS = ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include
LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a
- subroutine readpdb(lprint)
+ subroutine readpdb
C Read the PDB file and convert the peptide geometry into virtual-chain
C geometry.
- implicit none
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
- include 'COMMON.CONTROL'
+ include 'COMMON.FRAG'
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
- include 'COMMON.SBRIDGE'
- character*3 seq,atom,res
+ include 'COMMON.CONTROL'
+ integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
+ logical lprn /.false./,fail
+ double precision e1(3),e2(3),e3(3)
+ double precision dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
character*80 card
- double precision sccor(3,50)
- integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
- double precision dcj
- integer rescode,kkk,lll,icha,cou,kupa,iprzes
- logical lprint
+ double precision sccor(3,20)
+ integer rescode
+ efree_temp=0.0d0
ibeg=1
ishift1=0
+ ishift=0
+c write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
do
read (ipdbin,'(a80)',end=10) card
- if (card(:3).eq.'END') then
- goto 10
- else if (card(:3).eq.'TER') then
-C End current chain
-c ires_old=ires+1
- ires_old=ires+2
- itype(ires_old-1)=ntyp1
- itype(ires_old)=ntyp1
- ibeg=2
-c write (iout,*) "Chain ended",ires,ishift,ires_old
- call sccenter(ires,iii,sccor)
+c write (iout,'(a)') card
+ if (card(:5).eq.'HELIX') then
+ nhfrag=nhfrag+1
+ lsecondary=.true.
+ read(card(22:25),*) hfrag(1,nhfrag)
+ read(card(34:37),*) hfrag(2,nhfrag)
+ endif
+ if (card(:5).eq.'SHEET') then
+ nbfrag=nbfrag+1
+ lsecondary=.true.
+ read(card(24:26),*) bfrag(1,nbfrag)
+ read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc to be corrected !!!
+ bfrag(3,nbfrag)=bfrag(1,nbfrag)
+ bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
endif
+ if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+c Read free energy
+ if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
- read (card(14:16),'(a3)') atom
- if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(12:16),*) atom
+c write (iout,*) "! ",atom," !",ires
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+c write (iout,*) "ires",ires,ires-ishift+ishift1,
+c & " ires_old",ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
+c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
if (ibeg.eq.0) then
- call sccenter(ires,iii,sccor)
+c write (iout,*) "Calculating sidechain center iii",iii
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
endif
C Start new residue.
-c write (iout,'(a80)') card
- read (card(23:26),*) ires
- read (card(18:20),'(a3)') res
- if (ibeg.eq.1) then
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+c write (iout,*) "BEG ires",ires
ishift=ires-1
if (res.ne.'GLY' .and. res.ne. 'ACE') then
ishift=ishift-1
itype(1)=ntyp1
endif
-c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ires=ires-ishift+ishift1
+ ires_old=ires
+c write (iout,*) "ishift",ishift," ires",ires,
+c & " ires_old",ires_old
ibeg=0
- else if (ibeg.eq.2) then
-c Start a new chain
- ishift=-ires_old+ires-1
-c write (iout,*) "New chain started",ires,ishift
- ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
endif
- ires=ires-ishift
-c write (2,*) "ires",ires," ishift",ishift
- if (res.eq.'ACE') then
- ity=10
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires)=10
else
itype(ires)=rescode(ires,res,0)
endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+c write (iout,*) "ires_old",ires_old," ires",ires
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c ishift1=ishift1+1
+ endif
+c write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
- read(card(61:66),*) bfac(ires)
-c write (iout,'(2i3,2x,a,3f8.3,5x,f8.3)')
-c & ires,itype(ires),res,(c(j,ires),j=1,3),bfac(ires)
- iii=1
+c write (iout,*) "backbone ",atom
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)')
+ & ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
do j=1,3
sccor(j,iii)=c(j,ires)
enddo
- else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
- & atom(1:1).ne.'Q' .and. atom(1:2).ne.'1H' .and.
- & atom(1:2).ne.'2H' .and. atom(1:2).ne.'3H' .and.
- & atom.ne.'N ' .and. atom.ne.'C ' .and.
- & atom.ne.'OXT' ) then
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+c write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N' .and. atom.ne.'C' .and.
+ & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
+ & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c write (iout,*) "sidechain ",atom
iii=iii+1
-c write (iout,*) res,ires,iii,atom
read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
-c write (iout,'(3f8.3)') (sccor(j,iii),j=1,3)
endif
endif
enddo
- 10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
- nres=ires
- do i=2,nres-1
-c write (iout,*) i,itype(i)
-
- if (itype(i).eq.ntyp1) then
- if (itype(i+1).eq.ntyp1) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-C if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-C call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
-C if (fail) then
-C e2(1)=0.0d0
-C e2(2)=1.0d0
-C e2(3)=0.0d0
-C endif !fail
-C do j=1,3
-C c(j,i)=c(j,i-1)-1.9d0*e2(j)
-C enddo
-C else !unres_pdb
- do j=1,3
- dcj=(c(j,i-2)-c(j,i-3))/2.0
- c(j,i)=c(j,i-1)+dcj
- c(j,nres+i)=c(j,i)
- enddo
-C endif !unres_pdb
- else !itype(i+1).eq.ntyp1
-C if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
-C call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-C if (fail) then
-C e2(1)=0.0d0
-C e2(2)=1.0d0
-C e2(3)=0.0d0
-C endif
-C do j=1,3
-C c(j,i)=c(j,i+1)-1.9d0*e2(j)
-C enddo
-C else !unres_pdb
- do j=1,3
- dcj=(c(j,i+3)-c(j,i+2))/2.0
- c(j,i)=c(j,i+1)-dcj
- c(j,nres+i)=c(j,i)
- enddo
-C endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
- endif !itype.eq.ntyp1
- enddo
+ 10 continue
+#ifdef DEBUG
+ write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+ if (ires.eq.0) return
C Calculate the CM of the last side chain.
- call sccenter(ires,iii,sccor)
+ if (iii.gt.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+ nres=ires
nsup=nres
nstart_sup=1
if (itype(nres).ne.10) then
nres=nres+1
itype(nres)=ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+ enddo
+ else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
+ endif
endif
do i=2,nres-1
do j=1,3
if (itype(1).eq.ntyp1) then
nsup=nsup-1
nstart_sup=2
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-3.8d0*e2(j)
+ enddo
+ else
do j=1,3
- dcj=(c(j,4)-c(j,3))/2.0
+ dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
+ endif
endif
+C Copy the coordinates to reference coordinates
+c do i=1,2*nres
+c do j=1,3
+c cref(j,i)=c(j,i)
+c enddo
+c enddo
C Calculate internal coordinates.
- if (lprint) then
- write (iout,100)
+ if (lprn) then
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+C Calculate internal coordinates.
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
& ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
& (c(j,nres+ires),j=1,3)
- enddo
- endif
+ enddo
call int_from_cart(.true.,.false.)
- call flush(iout)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
do i=1,nres-1
do j=1,3
dc(j,i)=c(j,i+1)-c(j,i)
enddo
c call chainbuild
C Copy the coordinates to reference coordinates
- do i=1,nres
+ do i=1,2*nres
do j=1,3
cref(j,i)=c(j,i)
- cref(j,i+nres)=c(j,i+nres)
enddo
enddo
- 100 format ('Residue alpha-carbon coordinates ',
- & ' centroid coordinates'/
- 1 ' ', 6X,'X',7X,'Y',7X,'Z',
- & 12X,'X',7X,'Y',7X,'Z')
- 110 format (a,'(',i3,')',6f12.5)
+
+ do j=1,nbfrag
+ do i=1,4
+ bfrag(i,j)=bfrag(i,j)-ishift
+ enddo
+ enddo
+
+ do j=1,nhfrag
+ do i=1,2
+ hfrag(i,j)=hfrag(i,j)-ishift
+ enddo
+ enddo
ishift_pdb=ishift
return
end
c---------------------------------------------------------------------------
subroutine int_from_cart(lside,lprn)
- implicit none
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
- character*3 seq,atom,res
+ character*3 seq,res
+c character*5 atom
character*80 card
- double precision sccor(3,50)
+ dimension sccor(3,20)
integer rescode
- double precision dist,alpha,beta,di
- integer i,j,iti
logical lside,lprn
- if (lprn) then
+ if (lprn) then
write (iout,'(/a)')
& 'Internal coordinates calculated from crystal structure.'
if (lside) then
write (iout,'(8a)') ' Res ',' dvb',' Theta',
- & ' Phi',' Dsc_id',' Dsc',' Alpha',
- & ' Omega'
+ & ' Gamma',' Dsc_id',' Dsc',' Alpha',
+ & ' Beta '
else
write (iout,'(4a)') ' Res ',' dvb',' Theta',
- & ' Phi'
+ & ' Gamma'
endif
- endif
- do i=2,nres
+ endif
+ do i=1,nres-1
iti=itype(i)
-c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
- if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
- & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
+ if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
write (iout,'(a,i4)') 'Bad Cartesians for residue',i
- stop
+ctest stop
endif
- vbld(i)=dist(i-1,i)
- vbld_inv(i)=1.0d0/vbld(i)
- theta(i+1)=alpha(i-1,i,i+1)
+ vbld(i+1)=dist(i,i+1)
+ vbld_inv(i+1)=1.0d0/vbld(i+1)
+ if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
enddo
-c if (itype(1).eq.ntyp1) then
-c do j=1,3
-c c(j,1)=c(j,2)+(c(j,3)-c(j,4))
-c enddo
-c endif
-c if (itype(nres).eq.ntyp1) then
-c do j=1,3
-c c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
-c enddo
+c if (unres_pdb) then
+c if (itype(1).eq.ntyp1) then
+c theta(3)=90.0d0*deg2rad
+c phi(4)=180.0d0*deg2rad
+c vbld(2)=3.8d0
+c vbld_inv(2)=1.0d0/vbld(2)
+c endif
+c if (itype(nres).eq.ntyp1) then
+c theta(nres)=90.0d0*deg2rad
+c phi(nres)=180.0d0*deg2rad
+c vbld(nres)=3.8d0
+c vbld_inv(nres)=1.0d0/vbld(2)
+c endif
c endif
if (lside) then
do i=2,nres-1
do j=1,3
- c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
+ c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+ & +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
enddo
iti=itype(i)
di=dist(i,nres+i)
- vbld(i+nres)=di
+C 10/03/12 Adam: Correction for zero SC-SC bond length
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
+ & di=dsc(itype(i))
+ vbld(i+nres)=di
if (itype(i).ne.10) then
vbld_inv(i+nres)=1.0d0/di
else
alph(i)=alpha(nres+i,i,maxres2)
omeg(i)=beta(nres+i,i,maxres2,i+1)
endif
- if (iti.ne.10) then
- alph(i)=alpha(nres+i,i,maxres2)
- omeg(i)=beta(nres+i,i,maxres2,i+1)
- endif
- if (lprn)
- & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
- & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
- & rad2deg*alph(i),rad2deg*omeg(i)
+ if (lprn)
+ & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+ & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+ & rad2deg*alph(i),rad2deg*omeg(i)
enddo
else if (lprn) then
do i=2,nres
iti=itype(i)
write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
- & rad2deg*theta(i),rad2deg*phi(i)
+ & rad2deg*theta(i),rad2deg*phi(i)
enddo
endif
return
end
-c---------------------------------------------------------------------------
- subroutine sccenter(ires,nscat,sccor)
- implicit none
- include 'DIMENSIONS'
- include 'COMMON.CHAIN'
- integer ires,nscat,i,j
- double precision sccor(3,50),sccmj
- do j=1,3
- sccmj=0.0D0
- do i=1,nscat
- sccmj=sccmj+sccor(j,i)
- enddo
- dc(j,ires)=sccmj/nscat
- enddo
- return
- end
-c---------------------------------------------------------------------------
+c-------------------------------------------------------------------------------
subroutine sc_loc_geom(lprn)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
- include 'COMMON.SETUP'
double precision x_prime(3),y_prime(3),z_prime(3)
logical lprn
do i=1,nres-1
enddo
enddo
do i=2,nres-1
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i).ne.10) then
do j=1,3
dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
enddo
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
it=itype(i)
- if (it.ne.10 .and. itype(i).ne.ntyp1) then
+ if (it.ne.10) then
c
C Compute the axes of tghe local cartesian coordinates system; store in
c x_prime, y_prime and z_prime
x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
-c write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c write (iout,*) "y_prime",(y_prime(j),j=1,3)
call vecpr(x_prime,y_prime,z_prime)
-c write (iout,*) "z_prime",(z_prime(j),j=1,3)
c
C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
C to local coordinate system. Store in xx, yy, zz.
endif
enddo
if (lprn) then
- write (iout,*) "xxref,yyref,zzref"
do i=2,nres
iti=itype(i)
- write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
- & zzref(i)
+ write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
+ & yyref(i),zzref(i)
enddo
endif
return
end
c---------------------------------------------------------------------------
+ subroutine sccenter(ires,nscat,sccor)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ dimension sccor(3,20)
+ do j=1,3
+ sccmj=0.0D0
+ do i=1,nscat
+ sccmj=sccmj+sccor(j,i)
+ enddo
+ dc(j,ires)=sccmj/nscat
+ enddo
+ return
+ end
+c---------------------------------------------------------------------------
subroutine bond_regular
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
do i=1,nres-1
vbld(i+1)=vbl
vbld_inv(i+1)=1.0d0/vbld(i+1)
- vbld(i+1+nres)=dsc(iabs(itype(i+1)))
- vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
+ vbld(i+1+nres)=dsc(itype(i+1))
+ vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
c print *,vbld(i+1),vbld(i+1+nres)
enddo
return
end
-c---------------------------------------------------------------------------
subroutine readpdb_template(k)
-C Read the PDB file for read_constr_homology with read2sigma
+C Read the PDB file with gaps for read_constr_homology with read2sigma
C and convert the peptide geometry into virtual-chain geometry.
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
- include 'COMMON.SETUP'
integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
logical lprn /.false./,fail
double precision e1(3),e2(3),e3(3)
character*5 atom
character*80 card
double precision sccor(3,20)
- integer rescode,iterter(maxres)
- do i=1,maxres
- iterter(i)=0
- enddo
+ integer rescode
+ efree_temp=0.0d0
ibeg=1
ishift1=0
ishift=0
lsecondary=.false.
nhfrag=0
nbfrag=0
- do
+ do
read (ipdbin,'(a80)',end=10) card
- if (card(:3).eq.'END') then
- goto 10
- else if (card(:3).eq.'TER') then
-C End current chain
- ires_old=ires+2
- itype(ires_old-1)=ntyp1
- iterter(ires_old-1)=1
- itype(ires_old)=ntyp1
- iterter(ires_old)=1
- ibeg=2
-c write (iout,*) "Chain ended",ires,ishift,ires_old
- if (unres_pdb) then
- do j=1,3
- dc(j,ires)=sccor(j,iii)
- enddo
- else
- call sccenter(ires,iii,sccor)
- endif
- endif
+c write (iout,'(a)') card
+ if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
read (card(12:16),*) atom
c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
+c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
if (ibeg.eq.0) then
+c write (iout,*) "Calculating sidechain center iii",iii
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
ires_old=ires
c write (iout,*) "ishift",ishift," ires",ires,
c & " ires_old",ires_old
-c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
- ibeg=0
- else if (ibeg.eq.2) then
-c Start a new chain
- ishift=-ires_old+ires-1
- ires=ires_old+1
-c write (iout,*) "New chain started",ires,ishift
ibeg=0
else
ishift=ishift-(ires-ishift+ishift1-ires_old-1)
ires=ires-ishift+ishift1
endif
c write (iout,*) "ires_old",ires_old," ires",ires
-c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
c ishift1=ishift1+1
-c endif
+ endif
c write (2,*) "ires",ires," res ",res," ity",ity
if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
& res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+c write (iout,*) "backbone ",atom
#ifdef DEBUG
write (iout,'(2i3,2x,a,3f8.3)')
& ires,itype(ires),res,(c(j,ires),j=1,3)
endif
endif
enddo
- 10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
- nres=ires
- do i=2,nres-1
-c write (iout,*) i,itype(i),itype(i+1)
- if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
- if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif !fail
- do j=1,3
- c(j,i)=c(j,i-1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i-2)-c(j,i-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i-1)+dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- else !itype(i+1).eq.ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
- call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i+3)-c(j,i+2))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i+1)-dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
- endif !itype.eq.ntyp1
- enddo
+ 10 continue
+#ifdef DEBUG
+ write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+ if (ires.eq.0) return
C Calculate the CM of the last side chain.
+ if (iii.gt.0) then
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
else
call sccenter(ires,iii,sccor)
endif
+ endif
+ nres=ires
nsup=nres
nstart_sup=1
if (itype(nres).ne.10) then
nres=nres+1
itype(nres)=ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
- enddo
- else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
+ dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
endif
- endif
do i=2,nres-1
do j=1,3
c(j,i+nres)=dc(j,i)
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-1.9d0*e2(j)
+ c(j,1)=c(j,2)-3.8d0*e2(j)
enddo
else
do j=1,3
- dcj=(c(j,4)-c(j,3))/2.0
+ dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
endif
endif
-C Copy the coordinates to reference coordinates
-c do i=1,2*nres
-c do j=1,3
-c cref(j,i)=c(j,i)
-c enddo
-c enddo
C Calculate internal coordinates.
- if (out_template_coord) then
+ if (lprn) then
write (iout,'(/a)')
& "Cartesian coordinates of the reference structure"
write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
enddo
endif
C Calculate internal coordinates.
-c call int_from_cart1(.false.)
- call int_from_cart(.true.,.true.)
- call sc_loc_geom(.true.)
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
c & vbld_inv(i+nres)
enddo
- do i=1,nres
- do j=1,3
- cref(j,i)=c(j,i)
- cref(j,i+nres)=c(j,i+nres)
- enddo
- enddo
+c call chainbuild
+C Copy the coordinates to reference coordinates
do i=1,2*nres
do j=1,3
+ cref(j,i)=c(j,i)
chomo(j,i,k)=c(j,i)
enddo
enddo
+
+ ishift_pdb=ishift
return
end
-
-
cart2intgrad.o checkder_p.o contact_cp econstr_local.o econstr_qlike.o \
econstrq-PMF.o PMFprocess.o energy_p_new_barrier.o make_xx_list.o \
energy_p_new-sep_barrier.o gradient_p.o minimize_p.o sumsld.o \
- cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \
+ cored.o rmdd.o geomout.o readpdb.o int_from_cart.o regularize.o \
+ thread.o fitsq.o mcm.o \
mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o \
eigen.o blas.o add.o entmcm.o minim_mcmf.o \
together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \
${FC} ${FFLAGS} cinfo.f
${FC} ${OPT} ${object} cinfo.o ${LIBS} -o ${BIN}
+E0LL2Y_DFA: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
+ -DSPLITELE -DLANG0 -DFOURBODY -DDFA
+E0LL2Y_DFA: BIN = ~/bin/unres_ifort_MPICH-okeanos_E0LL2Y-HCD-DFA.exe
+E0LL2Y_DFA: ${object} dfa.o xdrf/libxdrf.a
+ gcc -o compinfo compinfo.c
+ ./compinfo | true
+ ${FC} ${FFLAGS} cinfo.f
+ ${FC} ${OPT} ${object} dfa.o cinfo.o ${LIBS} -o ${BIN}
+
NEWCORR: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DAMD64 -DUNRES -DISNAN -DMP -DMPI \
-DSPLITELE -DLANG0 -DNEWCORR -DCORRCD #-DFOURBODY #-DMYGAUSS #-DTIMING
NEWCORR: BIN = ~/bin/unres_ifort_MPICH-okeanos_SC-HCD.exe
readpdb.o : readpdb.F
${FC} ${FFLAGS2} ${CPPFLAGS} readpdb.F
+readpdb.o : readpdb.F
+ ${FC} ${FFLAGS2} ${CPPFLAGS} readpdb-mult.F
+
sumsld.o : sumsld.f
${FC} ${FFLAGS} ${CPPFLAGS} sumsld.f
--- /dev/null
+
+c------------------------------------------------------------------------
+ double precision function boxshift(x,boxsize)
+ implicit none
+ double precision x,boxsize
+ double precision xtemp
+ xtemp=dmod(x,boxsize)
+ if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+ boxshift=xtemp-boxsize
+ else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+ boxshift=xtemp+boxsize
+ else
+ boxshift=xtemp
+ endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine closest_img(xi,yi,zi,xj,yj,zj)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ integer xshift,yshift,zshift,subchap
+ double precision dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine to_box(xi,yi,zi)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ double precision xi,yi,zi
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0.0d0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0.0d0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0.0d0) zi=zi+boxzsize
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ double precision xi,yi,zi,sslipi,ssgradlipi
+ double precision fracinbuf
+ double precision sscalelip,sscagradlip
+
+ if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+ return
+ end
double precision app_(2,2),bpp_(2,2),rpp_(2,2)
integer ncont,icont(2,maxcont)
double precision econt(maxcont)
- integer xshift,yshift,zshift
+ double precision boxshift
*
* Load the constants of peptide bond - peptide bond interactions.
* Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
xmedi=xi+0.5*dxi
ymedi=yi+0.5*dyi
zmedi=zi+0.5*dzi
+ call to_box(xmedi,ymedi,zmedi)
c write (iout,*) "i",xmedi,ymedi,zmedi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
c write (iout,*) "i",xmedi,ymedi,zmedi
do 4 j=i+2,nct-1
c write (iout,*) "i",i," j",j
yj=c(2,j)+0.5*dyj
zj=c(3,j)+0.5*dzj
c write (iout,*) "j",xj,yj,zj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
-c write (iout,*) "j",xj,yj,zj
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-c write (iout,*) "dist",dsqrt(dist_init)
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
-c write (iout,*) "shift",xshift,yshift,zshift," dist_temp",
-c & dist_temp," dist_init",dist_init
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0/(xj*xj+yj*yj+zj*zj)
rmij=sqrt(rrmij)
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& sigij,r0ij,rcut,sss1,sssgrad1,sqrij
double precision sscale,sscagrad
+ double precision boxshift
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
sqrij=dsqrt(rrij)
eps0ij=eps(itypi,itypj)
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& sigij,r0ij,rcut,sqrij,sss1,sssgrad1
double precision sscale,sscagrad
+ double precision boxshift
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C Change 12/1/95
num_conti=0
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
C Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
sqrij=dsqrt(rij)
& fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
double precision sscale,sscagrad
+ double precision boxshift
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
& fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
double precision sscale,sscagrad
+ double precision boxshift
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
double precision sss1,sssgrad1
double precision sscale,sscagrad
+ double precision boxshift
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
integer itypi,itypj,itypi1,iint,ind,ikont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
double precision sscale,sscagrad
+ double precision boxshift
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
& xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision subchap,sss1,sssgrad1
+ double precision boxshift
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
xj=c(1,nres+j)
yj=c(2,nres+j)
zj=c(3,nres+j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- else if (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if (dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
- dxj=dc_norm(1,nres+j)
- dyj=dc_norm(2,nres+j)
- dzj=dc_norm(3,nres+j)
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- rij=dsqrt(rrij)
- sss1=sscale(1.0d0/rij,r_cut_int)
- if (sss1.eq.0.0d0) cycle
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
- if (sss.lt.1.0d0) then
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss1=sscale(1.0d0/rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+ if (sss.lt.1.0d0) then
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
- sssgrad=
+ sssgrad=
& sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
call sc_angular
include 'COMMON.CONTROL'
include "COMMON.SPLITELE"
logical lprn
- integer xshift,yshift,zshift
double precision evdw
integer itypi,itypj,itypi1,iint,ind,ikont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
- & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
- double precision subchap
+ double precision boxshift
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
xj=c(1,nres+j)
yj=c(2,nres+j)
zj=c(3,nres+j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
& +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
& +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
& sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
& xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision boxshift
evdw=0.0D0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
+ call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
c dsci_inv=dsc_inv(itypi)
dsci_inv=vbld_inv(i+nres)
C
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
num_conti=0
call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
#ifdef FOURBODY
num_conti=num_cont_hb(i)
#endif
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
#ifdef FOURBODY
num_conti=num_cont_hb(i)
double precision sss1,sssgrad1
double precision sscale,sscagrad
double precision scalar
+ double precision boxshift
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
-
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0D0/rij
rij=dsqrt(rij)
double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
& dist_temp, dist_init,sss_grad
double precision sscale,sscagrad
+ double precision boxshift
integer ikont
evdw1=0.0D0
C print *,"WCHODZE"
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
c & ' ielend',ielend_vdw(i)
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0D0/rij
rij=dsqrt(rij)
logical lprint_short
common /shortcheck/ lprint_short
double precision ggg(3)
- integer xshift,yshift,zshift
integer i,iint,j,k,iteli,itypj,subchap
double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
& fac,e1,e2,rij
double precision evdw2,evdw2_14,evdwij
- double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
- & dist_temp, dist_init
double precision sscale,sscagrad
+ double precision boxshift
integer ikont
if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
evdw2=0.0D0
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
c do iint=1,nscp_gr(i)
yj=c(2,j)
zj=c(3,j)
c corrected by AL
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
-c end correction
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
& dist_temp, dist_init
double precision ggg(3)
double precision sscale,sscagrad
+ double precision boxshift
evdw2=0.0D0
evdw2_14=0.0d0
cd print '(a)','Enter ESCP'
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
c if (lprint_short)
c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
yj=c(2,j)
zj=c(3,j)
c corrected by AL
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
-c end correction
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-c if (lprint_short) then
-c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
-c write (iout,*) "dist_init",dsqrt(dist_init)
-c endif
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
-c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
if (nfgtasks.gt.1) then
call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
endif
+c write (iout,*) "itime_mat",itime_mat," imatupdate",imatupdate
if (mod(itime_mat,imatupdate).eq.0) then
call make_SCp_inter_list
call make_SCSC_inter_list
& sigij,r0ij,rcut,sqrij,sss1,sssgrad1
double precision fcont,fprimcont
double precision sscale,sscagrad
+ double precision boxshift
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C Change 12/1/95
num_conti=0
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
C Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
rrij=1.0D0/rij
& fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
double precision sscale,sscagrad
+ double precision boxshift
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi,
& sss1,sssgrad1
double precision sscale,sscagrad
+ double precision boxshift
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
include 'COMMON.SPLITELE'
include 'COMMON.SBRIDGE'
logical lprn
- integer xshift,yshift,zshift,subchap
double precision evdw
integer itypi,itypj,itypi1,iint,ind,ikont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
- & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision boxshift
evdw=0.0D0
ccccc energy_dec=.false.
C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
-C Return atom into box, boxxsize is size of box in x dimension
-c 134 continue
-c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
-c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
-c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
-c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
-c go to 134
-c endif
-c 135 continue
-c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
-c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
-c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
-c & (yi.lt.((yshift-0.5d0)*boxysize))) then
-c go to 135
-c endif
-c 136 continue
-c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
-c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
-C Condition for being inside the proper box
-c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
-c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
-c go to 136
-c endif
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
C define scaling factor for lipids
C if (positi.le.0) positi=positi+boxzsize
C print *,i
C first for peptide groups
c for each residue check if it is in lipid or lipid water border area
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
-
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
C xi=xi+xshift*boxxsize
C yi=yi+yshift*boxysize
C zi=zi+zshift*boxzsize
if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
& 'evdw',i,j,evdwij,' ss'
C triple bond artifac removal
- do k=j+1,iend(i,iint)
+ do k=j+1,iend(i,iint)
C search over all next residues
- if (dyn_ss_mask(k)) then
+ if (dyn_ss_mask(k)) then
C check if they are cysteins
C write(iout,*) 'k=',k
c write(iout,*) "PRZED TRI", evdwij
- evdwij_przed_tri=evdwij
- call triple_ssbond_ene(i,j,k,evdwij)
+ evdwij_przed_tri=evdwij
+ call triple_ssbond_ene(i,j,k,evdwij)
c if(evdwij_przed_tri.ne.evdwij) then
c write (iout,*) "TRI:", evdwij, evdwij_przed_tri
c endif
c write(iout,*) "PO TRI", evdwij
C call the energy function that removes the artifical triple disulfide
C bond the soubroutine is located in ssMD.F
- evdw=evdw+evdwij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ evdw=evdw+evdwij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
& 'evdw',i,j,evdwij,'tss'
- endif!dyn_ss_mask(k)
- enddo! k
+ endif!dyn_ss_mask(k)
+ enddo! k
ELSE
- ind=ind+1
- itypj=iabs(itype(j))
- if (itypj.eq.ntyp1) cycle
+ ind=ind+1
+ itypj=iabs(itype(j))
+ if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
- dscj_inv=vbld_inv(j+nres)
+ dscj_inv=vbld_inv(j+nres)
c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
c & 1.0d0/vbld(j+nres)
c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
- sig0ij=sigma(itypi,itypj)
- chi1=chi(itypi,itypj)
- chi2=chi(itypj,itypi)
- chi12=chi1*chi2
- chip1=chip(itypi)
- chip2=chip(itypj)
- chip12=chip1*chip2
- alf1=alp(itypi)
- alf2=alp(itypj)
- alf12=0.5D0*(alf1+alf2)
+ sig0ij=sigma(itypi,itypj)
+ chi1=chi(itypi,itypj)
+ chi2=chi(itypj,itypi)
+ chi12=chi1*chi2
+ chip1=chip(itypi)
+ chip2=chip(itypj)
+ chip12=chip1*chip2
+ alf1=alp(itypi)
+ alf2=alp(itypj)
+ alf12=0.5D0*(alf1+alf2)
C For diagnostics only!!!
c chi1=0.0D0
c chi2=0.0D0
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)
- yj=c(2,nres+j)
- zj=c(3,nres+j)
-C Return atom J into box the original box
-c 137 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 137
-c endif
-c 138 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 138
-c endif
-c 139 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 139
-c endif
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot)
- &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zj.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
C write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)
C if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
C if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
C print *,sslipi,sslipj,bordlipbot,zi,zj
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
- dxj=dc_norm(1,nres+j)
- dyj=dc_norm(2,nres+j)
- dzj=dc_norm(3,nres+j)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
C xj=xj-xi
C yj=yj-yi
C zj=zj-zi
c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
c write (iout,*) "j",j," dc_norm",
c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- rij=dsqrt(rrij)
- sss=sscale(1.0d0/rij,r_cut_int)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/rij,r_cut_int)
c write (iout,'(a7,4f8.3)')
c & "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
- if (sss.eq.0.0d0) cycle
- sssgrad=sscagrad(1.0d0/rij,r_cut_int)
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
- call sc_angular
- sigsq=1.0D0/sigsq
- sig=sig0ij*dsqrt(sigsq)
- rij_shift=1.0D0/rij-sig+sig0ij
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+sig0ij
c for diagnostics; uncomment
c rij_shift=1.2*sig0ij
C I hate to put IF's in the loops, but here don't have another choice!!!!
- if (rij_shift.le.0.0D0) then
- evdw=1.0D20
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
- return
- endif
- sigder=-sig*sigsq
+ return
+ endif
+ sigder=-sig*sigsq
c---------------------------------------------------------------
- rij_shift=1.0D0/rij_shift
- fac=rij_shift**expon
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
C here to start with
C if (c(i,3).gt.
- faclip=fac
- e1=fac*fac*aa
- e2=fac*bb
- evdwij=eps1*eps2rt*eps3rt*(e1+e2)
- eps2der=evdwij*eps3rt
- eps3der=evdwij*eps2rt
+ faclip=fac
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
C &((sslipi+sslipj)/2.0d0+
C &(2.0d0-sslipi-sslipj)/2.0d0)
c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
- evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
- if (lprn) then
- sigm=dabs(aa/bb)**(1.0D0/6.0D0)
- epsi=bb**2/aa
- write (iout,'(2(a3,i3,2x),17(0pf7.3))')
- & restyp(itypi),i,restyp(itypj),j,
- & epsi,sigm,chi1,chi2,chip1,chip2,
- & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
- & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
- & evdwij
- endif
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij*sss
+ if (lprn) then
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ & restyp(itypi),i,restyp(itypj),j,
+ & epsi,sigm,chi1,chi2,chip1,chip2,
+ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+ & evdwij
+ endif
- if (energy_dec) write (iout,'(a,2i5,3f10.5)')
+ if (energy_dec) write (iout,'(a,2i5,3f10.5)')
& 'r sss evdw',i,j,1.0d0/rij,sss,evdwij
C Calculate gradient components.
- e1=e1*eps1*eps2rt**2*eps3rt**2
- fac=-expon*(e1+evdwij)*rij_shift
- sigder=fac*sigder
- fac=rij*fac
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac
c print '(2i4,6f8.4)',i,j,sss,sssgrad*
c & evdwij,fac,sigma(itypi,itypj),expon
- fac=fac+evdwij*sssgrad/sss*rij
+ fac=fac+evdwij*sssgrad/sss*rij
c fac=0.0d0
C Calculate the radial part of the gradient
- gg_lipi(3)=eps1*(eps2rt*eps2rt)
- & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
- & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
- & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
- gg_lipj(3)=ssgradlipj*gg_lipi(3)
- gg_lipi(3)=gg_lipi(3)*ssgradlipi
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
C gg_lipi(3)=0.0d0
C gg_lipj(3)=0.0d0
- gg(1)=xj*fac
- gg(2)=yj*fac
- gg(3)=zj*fac
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
C Calculate angular part of the gradient.
c call sc_grad_scale(sss)
- call sc_grad
+ call sc_grad
ENDIF ! dyn_ss
c enddo ! j
c enddo ! iint
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.SPLITELE'
- integer xshift,yshift,zshift,subchap
+ double precision boxshift
integer icall
common /srutu/ icall
logical lprn
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
& xi,yi,zi,fac_augm,e_augm
double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
- & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip,sssgrad1
+ & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip,sssgrad1
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
evdw=0.0D0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
C define scaling factor for lipids
C if (positi.le.0) positi=positi+boxzsize
C print *,i
C first for peptide groups
c for each residue check if it is in lipid or lipid water border area
- if ((zi.gt.bordlipbot)
- &.and.(zi.lt.bordliptop)) then
-C the energy transfer exist
- if (zi.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
- sslipi=sscalelip(fracinbuf)
- ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipi=1.0d0
- ssgradlipi=0.0
- endif
- else
- sslipi=0.0d0
- ssgradlipi=0.0
- endif
-
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
-C xj=c(1,nres+j)-xi
-C yj=c(2,nres+j)-yi
-C zj=c(3,nres+j)-zi
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.gt.bordlipbot)
- &.and.(zj.lt.bordliptop)) then
-C the energy transfer exist
- if (zj.lt.buflipbot) then
-C what fraction I am in
- fracinbuf=1.0d0-
- & ((zj-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zj.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
- else
- sslipj=1.0d0
- ssgradlipj=0.0
- endif
- else
- sslipj=0.0d0
- ssgradlipj=0.0
- endif
- aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
- bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
C if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5')
C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
C write(iout,*) "tu,", i,j,aa,bb,aa_lip(itypi,itypj),sslipi,sslipj
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
- dxj=dc_norm(1,nres+j)
- dyj=dc_norm(2,nres+j)
- dzj=dc_norm(3,nres+j)
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- rij=dsqrt(rrij)
- sss=sscale(1.0d0/rij,r_cut_int)
- if (sss.eq.0.0d0) cycle
- sssgrad=sscagrad(1.0d0/rij,r_cut_int)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+ rij=dsqrt(rrij)
+ sss=sscale(1.0d0/rij,r_cut_int)
+ if (sss.eq.0.0d0) cycle
+ sssgrad=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
- call sc_angular
- sigsq=1.0D0/sigsq
- sig=sig0ij*dsqrt(sigsq)
- rij_shift=1.0D0/rij-sig+r0ij
+ call sc_angular
+ sigsq=1.0D0/sigsq
+ sig=sig0ij*dsqrt(sigsq)
+ rij_shift=1.0D0/rij-sig+r0ij
C I hate to put IF's in the loops, but here don't have another choice!!!!
- if (rij_shift.le.0.0D0) then
- evdw=1.0D20
- return
- endif
- sigder=-sig*sigsq
+ if (rij_shift.le.0.0D0) then
+ evdw=1.0D20
+ return
+ endif
+ sigder=-sig*sigsq
c---------------------------------------------------------------
- rij_shift=1.0D0/rij_shift
- fac=rij_shift**expon
- e1=fac*fac*aa
- e2=fac*bb
- evdwij=eps1*eps2rt*eps3rt*(e1+e2)
- eps2der=evdwij*eps3rt
- eps3der=evdwij*eps2rt
- fac_augm=rrij**expon
- e_augm=augm(itypi,itypj)*fac_augm
- evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij+e_augm
- if (lprn) then
- sigm=dabs(aa/bb)**(1.0D0/6.0D0)
- epsi=bb**2/aa
- write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+ rij_shift=1.0D0/rij_shift
+ fac=rij_shift**expon
+ e1=fac*fac*aa
+ e2=fac*bb
+ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+ eps2der=evdwij*eps3rt
+ eps3der=evdwij*eps2rt
+ fac_augm=rrij**expon
+ e_augm=augm(itypi,itypj)*fac_augm
+ evdwij=evdwij*eps2rt*eps3rt
+ evdw=evdw+evdwij+e_augm
+ if (lprn) then
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
+ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
& chi1,chi2,chip1,chip2,
& eps1,eps2rt**2,eps3rt**2,
& om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
& evdwij+e_augm
- endif
+ endif
C Calculate gradient components.
- e1=e1*eps1*eps2rt**2*eps3rt**2
- fac=-expon*(e1+evdwij)*rij_shift
- sigder=fac*sigder
- fac=rij*fac-2*expon*rrij*e_augm
- fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ fac=-expon*(e1+evdwij)*rij_shift
+ sigder=fac*sigder
+ fac=rij*fac-2*expon*rrij*e_augm
+ fac=fac+(evdwij+e_augm)*sssgrad/sss*rij
C Calculate the radial part of the gradient
- gg(1)=xj*fac
- gg(2)=yj*fac
- gg(3)=zj*fac
+ gg(1)=xj*fac
+ gg(2)=yj*fac
+ gg(3)=zj*fac
C Calculate angular part of the gradient.
c call sc_grad_scale(sss)
- call sc_grad
+ call sc_grad
c enddo ! j
c enddo ! iint
enddo ! i
include 'COMMON.IOUNITS'
c include 'COMMON.CONTACTS'
dimension gg(3)
+ double precision boxshift
cd print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
evdw=0.0D0
c do i=iatsc_s,iatsc_e
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ call to_box(xi,yi,zi)
C
C Calculate SC interaction energy.
C
c do j=istart(i,iint),iend(i,iint)
itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ xj=boxshift(c(1,nres+j)-xi,boxxsize)
+ yj=boxshift(c(2,nres+j)-yi,boxysize)
+ zj=boxshift(c(3,nres+j)-zi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
r0ij=r0(itypi,itypj)
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
dimension ggg(3)
- integer xshift,yshift,zshift
+ double precision boxshift
C write(iout,*) 'In EELEC_soft_sphere'
ees=0.0D0
evdw1=0.0D0
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
do j=ielstart(i),ielend(i)
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
rij=xj*xj+yj*yj+zj*zj
sss=sscale(sqrt(rij),r_cut_int)
sssgrad=sscagrad(sqrt(rij),r_cut_int)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
c & (zmedi.lt.((-0.5d0)*boxzsize))) then
c go to 196
c endif
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
-
+ call to_box(xmedi,ymedi,zmedi)
#ifdef FOURBODY
num_conti=num_cont_hb(i)
#endif
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0) zmedi=zmedi+boxzsize
+ call to_box(xmedi,ymedi,zmedi)
C xmedi=xmedi+xshift*boxxsize
C ymedi=ymedi+yshift*boxysize
C zmedi=zmedi+zshift*boxzsize
C do j=16,17
C write (iout,*) i,j
C if (j.le.1) cycle
- if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
+ if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
C changes suggested by Ana to avoid out of bounds
c & .or.((j+2).gt.nres)
c & .or.((j-1).le.0)
c & .or.itype(j+2).eq.ntyp1
c & .or.itype(j-1).eq.ntyp1
&) cycle
- call eelecij(i,j,ees,evdw1,eel_loc)
+ call eelecij(i,j,ees,evdw1,eel_loc)
c enddo ! j
#ifdef FOURBODY
num_cont_hb(i)=num_conti
& ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
& ecosgp,ecosam,ecosbm,ecosgm,ghalf,rlocshield
double precision a22,a23,a32,a33,geel_loc_ij,geel_loc_ji
- double precision dist_init,xj_safe,yj_safe,zj_safe,
- & xj_temp,yj_temp,zj_temp,dist_temp,xmedi,ymedi,zmedi
+ double precision xmedi,ymedi,zmedi
double precision sscale,sscagrad,scalar
-
+ double precision boxshift
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /1.0d0/
double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
& 0.0d0,1.0d0,0.0d0,
& 0.0d0,0.0d0,1.0d0/
- integer xshift,yshift,zshift
c time00=MPI_Wtime()
cd write (iout,*) "eelecij",i,j
c ind=ind+1
xj=c(1,j)+0.5D0*dxj
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xmedi,boxxsize)
+ yj=boxshift(yj-ymedi,boxysize)
+ zj=boxshift(zj-zmedi,boxzsize)
C if ((i+3).lt.j) then !this condition keeps for turn3 and turn4 not subject to PBC
c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-C Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
-c endif
-C endif !endPBC condintion
-C xj=xj-xmedi
-C yj=yj-ymedi
-C zj=zj-zmedi
rij=xj*xj+yj*yj+zj*zj
sss=sscale(dsqrt(rij),r_cut_int)
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
dimension ggg(3)
- integer xshift,yshift,zshift
+ double precision boxshift
evdw2=0.0D0
evdw2_14=0.0d0
r0_scp=4.5d0
c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
c go to 136
c endif
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ call to_box(xi,yi,zi)
C xi=xi+xshift*boxxsize
C yi=yi+yshift*boxysize
C zi=zi+zshift*boxzsize
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
-c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-cC Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
-c c endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
C xj=xj-xi
C yj=yj-yi
C zj=zj-zi
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
-cgrad if (j.lt.i) then
-cd write (iout,*) 'j<i'
-C Uncomment following three lines for SC-p interactions
-c do k=1,3
-c gradx_scp(k,j)=gradx_scp(k,j)+ggg(k)
-c enddo
-cgrad else
-cd write (iout,*) 'j>i'
-cgrad do k=1,3
-cgrad ggg(k)=-ggg(k)
-C Uncomment following line for SC-p interactions
-c gradx_scp(k,j)=gradx_scp(k,j)-ggg(k)
-cgrad enddo
-cgrad endif
-cgrad do k=1,3
-cgrad gvdwc_scp(k,i)=gvdwc_scp(k,i)-0.5D0*ggg(k)
-cgrad enddo
-cgrad kstart=min0(i+1,j)
-cgrad kend=max0(i-1,j-1)
-cd write (iout,*) 'i=',i,' j=',j,' kstart=',kstart,' kend=',kend
-cd write (iout,*) ggg(1),ggg(2),ggg(3)
-cgrad do k=kstart,kend
-cgrad do l=1,3
-cgrad gvdwc_scp(l,k)=gvdwc_scp(l,k)-ggg(l)
-cgrad enddo
-cgrad enddo
do k=1,3
gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
- integer xshift,yshift,zshift
double precision ggg(3)
integer i,iint,j,k,iteli,itypj,subchap,ikont
double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
& fac,e1,e2,rij
double precision evdw2,evdw2_14,evdwij
- double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
- & dist_temp, dist_init
double precision sscale,sscagrad
+ double precision boxshift
evdw2=0.0D0
evdw2_14=0.0d0
c print *,boxxsize,boxysize,boxzsize,'wymiary pudla'
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
-c xi=xi+xshift*boxxsize
-c yi=yi+yshift*boxysize
-c zi=zi+zshift*boxzsize
-c print *,xi,yi,zi,'polozenie i'
-C Return atom into box, boxxsize is size of box in x dimension
-c 134 continue
-c if (xi.gt.((xshift+0.5d0)*boxxsize)) xi=xi-boxxsize
-c if (xi.lt.((xshift-0.5d0)*boxxsize)) xi=xi+boxxsize
-C Condition for being inside the proper box
-c if ((xi.gt.((xshift+0.5d0)*boxxsize)).or.
-c & (xi.lt.((xshift-0.5d0)*boxxsize))) then
-c go to 134
-c endif
-c 135 continue
-c print *,xi,boxxsize,"pierwszy"
-
-c if (yi.gt.((yshift+0.5d0)*boxysize)) yi=yi-boxysize
-c if (yi.lt.((yshift-0.5d0)*boxysize)) yi=yi+boxysize
-C Condition for being inside the proper box
-c if ((yi.gt.((yshift+0.5d0)*boxysize)).or.
-c & (yi.lt.((yshift-0.5d0)*boxysize))) then
-c go to 135
-c endif
-c 136 continue
-c if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
-c if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
-C Condition for being inside the proper box
-c if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
-c & (zi.lt.((zshift-0.5d0)*boxzsize))) then
-c go to 136
-c endif
+ call to_box(xi,yi,zi)
c do iint=1,nscp_gr(i)
c do j=iscpstart(i,iint),iscpend(i,iint)
xj=c(1,j)
yj=c(2,j)
zj=c(3,j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
-c 174 continue
-c if (xj.gt.((0.5d0)*boxxsize)) xj=xj-boxxsize
-c if (xj.lt.((-0.5d0)*boxxsize)) xj=xj+boxxsize
-C Condition for being inside the proper box
-c if ((xj.gt.((0.5d0)*boxxsize)).or.
-c & (xj.lt.((-0.5d0)*boxxsize))) then
-c go to 174
-c endif
-c 175 continue
-c if (yj.gt.((0.5d0)*boxysize)) yj=yj-boxysize
-c if (yj.lt.((-0.5d0)*boxysize)) yj=yj+boxysize
-cC Condition for being inside the proper box
-c if ((yj.gt.((0.5d0)*boxysize)).or.
-c & (yj.lt.((-0.5d0)*boxysize))) then
-c go to 175
-c endif
-c 176 continue
-c if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
-c if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
-C Condition for being inside the proper box
-c if ((zj.gt.((0.5d0)*boxzsize)).or.
-c & (zj.lt.((-0.5d0)*boxzsize))) then
-c go to 176
-c endif
-CHERE IS THE CALCULATION WHICH MIRROR IMAGE IS THE CLOSEST ONE
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ call to_box(xj,yj,zj)
+ xj=boxshift(xj-xi,boxxsize)
+ yj=boxshift(yj-yi,boxysize)
+ zj=boxshift(zj-zi,boxzsize)
c print *,xj,yj,zj,'polozenie j'
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
c print *,rrij
endif
return
end
+c------------------------------------------------------------------------
+ double precision function boxshift(x,boxsize)
+ implicit none
+ double precision x,boxsize
+ double precision xtemp
+ xtemp=dmod(x,boxsize)
+ if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
+ boxshift=xtemp-boxsize
+ else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
+ boxshift=xtemp+boxsize
+ else
+ boxshift=xtemp
+ endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine closest_img(xi,yi,zi,xj,yj,zj)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ integer xshift,yshift,zshift,subchap
+ double precision dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine to_box(xi,yi,zi)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ double precision xi,yi,zi
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0.0d0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0.0d0) yi=yi+boxysize
+ zi=dmod(zi,boxzsize)
+ if (zi.lt.0.0d0) zi=zi+boxzsize
+ return
+ end
+c--------------------------------------------------------------------------
+ subroutine lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ double precision xi,yi,zi,sslipi,ssgradlipi
+ double precision fracinbuf
+ double precision sscalelip,sscagradlip
+
+ if ((zi.gt.bordlipbot).and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+ return
+ end
integer ioverlap(maxres),ioverlap_last
data redfac /0.5D0/
- write (iout,*) "overlap_sc_list"
+c write (iout,*) "overlap_sc_list"
ioverlap_last=0
C Check for SC-SC overlaps and mark residues
c print *,'>>overlap_sc nnt=',nnt,' nct=',nct
ind=0
+c write(iout,*) "iatsc_s",iatsc_s," iatsc_e",iatsc_e
do i=iatsc_s,iatsc_e
itypi=iabs(itype(i))
itypi1=iabs(itype(i+1))
else
rcomp=sigma(itypi,itypj)
endif
-c print '(2(a3,2i3),a3,2f10.5)',
+c write (iout,'(2(a3,2i3),a3,2f10.5)'),
c & ' i=',i,iti,' j=',j,itj,' d=',dist(nres+i,nres+j)
c & ,rcomp
xj=c(1,nres+j)-xi
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
call sc_angular
+c write (iout,*) "dxj",dxj," dyj",dyj," dzj",dzj
+c write (iout,*) "erij",erij
+c write (iout,*) "om1",om1," om2",om2," om12",om12,
+c & " faceps1",faceps1," eps1",eps1
+c write (iout,*) "sigsq",sigsq
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
rij_shift=1.0D0/rij-sig+sig0ij
-
-ct if ( 1.0/rij .lt. redfac*rcomp .or.
-ct & rij_shift.le.0.0D0 ) then
-c write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
+c write (iout,*) "rij_shift",rij_shift
+c if ( 1.0/rij .lt. redfac*rcomp .or.
+c & rij_shift.le.0.0D0 ) then
+c write (iout,'(a,i3,a,i3,a,f10.5,a,3f10.5)')
c & 'overlap SC-SC: i=',i,' j=',j,
c & ' dist=',dist(nres+i,nres+j),' rcomp=',
c & rcomp,1.0/rij,rij_shift
if ( rij_shift.le.0.0D0 ) then
+c write (iout,*) "overlap",i,j
ioverlap_last=ioverlap_last+1
ioverlap(ioverlap_last)=i
do k=1,ioverlap_last-1
--- /dev/null
+ subroutine int_from_cart(lside,lprn)
+ implicit none
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.NAMES'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SETUP'
+ double precision dist,alpha,beta
+ character*3 seq,atom,res
+ character*80 card
+ double precision sccor(3,50)
+ integer rescode
+ logical lside,lprn
+ integer i,j,iti
+ double precision di,cosfac2,sinfac2,cosfac,sinfac
+#ifdef MPI
+ if(me.eq.king.or..not.out1file)then
+#endif
+ if (lprn) then
+ write (iout,'(/a)')
+ & 'Internal coordinates calculated from crystal structure.'
+ if (lside) then
+ write (iout,'(8a)') ' Res ',' dvb',' Theta',
+ & ' Phi',' Dsc_id',' Dsc',' Alpha',
+ & ' Omega'
+ else
+ write (iout,'(4a)') ' Res ',' dvb',' Theta',
+ & ' Phi'
+ endif
+ endif
+#ifdef MPI
+ endif
+#endif
+ do i=1,nres-1
+ iti=itype(i)
+ if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
+ & (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
+ write (iout,'(a,i4)') 'Bad Cartesians for residue',i
+ctest stop
+ endif
+ vbld(i+1)=dist(i,i+1)
+ vbld_inv(i+1)=1.0d0/vbld(i+1)
+c write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
+c & vbld_inv(i+1)
+ if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
+ if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+ enddo
+c if (unres_pdb) then
+c if (itype(1).eq.21) then
+c theta(3)=90.0d0*deg2rad
+c phi(4)=180.0d0*deg2rad
+c vbld(2)=3.8d0
+c vbld_inv(2)=1.0d0/vbld(2)
+c endif
+c if (itype(nres).eq.21) then
+c theta(nres)=90.0d0*deg2rad
+c phi(nres)=180.0d0*deg2rad
+c vbld(nres)=3.8d0
+c vbld_inv(nres)=1.0d0/vbld(2)
+c endif
+c endif
+c print *,"A TU2"
+ if (lside) then
+ do i=2,nres-1
+ do j=1,3
+ c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+ & +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
+ enddo
+ iti=itype(i)
+ di=dist(i,nres+i)
+ vbld(i+nres)=di
+ if (itype(i).ne.10) then
+ vbld_inv(i+nres)=1.0d0/di
+ else
+ vbld_inv(i+nres)=0.0d0
+ endif
+ if (iti.ne.10) then
+ alph(i)=alpha(nres+i,i,maxres2)
+ omeg(i)=beta(nres+i,i,maxres2,i+1)
+ endif
+ if(me.eq.king.or..not.out1file)then
+ if (lprn)
+ & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+ & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+ & rad2deg*alph(i),rad2deg*omeg(i)
+ endif
+ enddo
+ if (lprn) then
+ i=nres
+ iti=itype(nres)
+ write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+ & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+ & rad2deg*alph(i),rad2deg*omeg(i)
+ endif
+ else if (lprn) then
+ do i=2,nres
+ iti=itype(i)
+ if(me.eq.king.or..not.out1file)
+ & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
+ & rad2deg*theta(i),rad2deg*phi(i)
+ enddo
+ endif
+ return
+ end
+c-------------------------------------------------------------------------------
+ subroutine sc_loc_geom(lprn)
+ implicit none
+ include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+#endif
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.NAMES'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SETUP'
+ double precision x_prime(3),y_prime(3),z_prime(3)
+ logical lprn
+ integer i,j,it
+ double precision xx,yy,zz,cosfac,cosfac2,sinfac,sinfac2
+ do i=1,nres-1
+ do j=1,3
+ dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
+ enddo
+c write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
+c & " vbld",vbld_inv(i+1)
+ enddo
+ do i=2,nres-1
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
+ enddo
+c write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
+c & " vbld",vbld_inv(i+nres)
+ else
+ do j=1,3
+ dc_norm(j,i+nres)=0.0d0
+ enddo
+ endif
+ enddo
+ do i=2,nres-1
+ costtab(i+1) =dcos(theta(i+1))
+ sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
+ cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
+ sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
+ cosfac2=0.5d0/(1.0d0+costtab(i+1))
+ cosfac=dsqrt(cosfac2)
+ sinfac2=0.5d0/(1.0d0-costtab(i+1))
+ sinfac=dsqrt(sinfac2)
+ it=itype(i)
+c write (iout,*) "i",i," costab",costtab(i+1),
+c & " sintab",sinttab(i+1)
+c write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
+c write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
+ if (it.ne.10 .and. itype(i).ne.ntyp1) then
+c
+C Compute the axes of tghe local cartesian coordinates system; store in
+c x_prime, y_prime and z_prime
+c
+ do j=1,3
+ x_prime(j) = 0.00
+ y_prime(j) = 0.00
+ z_prime(j) = 0.00
+ enddo
+ do j = 1,3
+ x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
+ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
+ enddo
+c write (iout,*) "x_prime",(x_prime(j),j=1,3)
+c write (iout,*) "y_prime",(y_prime(j),j=1,3)
+ call vecpr(x_prime,y_prime,z_prime)
+c write (iout,*) "z_prime",(z_prime(j),j=1,3)
+c
+C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
+C to local coordinate system. Store in xx, yy, zz.
+c
+ xx=0.0d0
+ yy=0.0d0
+ zz=0.0d0
+ do j = 1,3
+ xx = xx + x_prime(j)*dc_norm(j,i+nres)
+ yy = yy + y_prime(j)*dc_norm(j,i+nres)
+ zz = zz + z_prime(j)*dc_norm(j,i+nres)
+ enddo
+
+ xxref(i)=xx
+ yyref(i)=yy
+ zzref(i)=zz
+ else
+ xxref(i)=0.0d0
+ yyref(i)=0.0d0
+ zzref(i)=0.0d0
+ endif
+ enddo
+ if (lprn) then
+#ifdef MPI
+ if (me.eq.king.or..not.out1file) then
+#endif
+ write (iout,*) "xxref,yyref,zzref"
+ do i=2,nres
+ write (iout,'(a3,i4,3f10.5)')
+ & restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
+ enddo
+#ifdef MPI
+ endif
+#endif
+ endif
+ return
+ end
+c---------------------------------------------------------------------------
+ subroutine sccenter(ires,nscat,sccor)
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ integer i,j,ires,nscat
+ double precision sccor(3,50)
+ double precision sccmj
+ do j=1,3
+ sccmj=0.0D0
+ do i=1,nscat
+ sccmj=sccmj+sccor(j,i)
+ enddo
+ dc(j,ires)=sccmj/nscat
+ enddo
+ return
+ end
+c---------------------------------------------------------------------------
+ subroutine bond_regular
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.INTERACT'
+ include 'COMMON.CHAIN'
+ integer i,i1,i2
+ do i=1,nres-1
+ vbld(i+1)=vbl
+ vbld_inv(i+1)=vblinv
+ vbld(i+1+nres)=dsc(iabs(itype(i+1)))
+ vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
+c print *,vbld(i+1),vbld(i+1+nres)
+ enddo
+c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
+ do i=1,nchain
+ i1=chain_border(1,i)
+ i2=chain_border(2,i)
+ if (i1.gt.1) then
+ vbld(i1)=vbld(i1)/2
+ vbld_inv(i1)=vbld_inv(i1)*2
+ endif
+ if (i2.lt.nres) then
+ vbld(i2+1)=vbld(i2+1)/2
+ vbld_inv(i2+1)=vbld_inv(i2+1)*2
+ endif
+ enddo
+ return
+ end
endif
enddo
nvarx=k
- write (iout,*) "Variables set up nvarx",nvarx
- write (iout,*) "Before energy minimization"
+ call chainbuild_cart
call etotal(energia(0))
call enerprint(energia(0))
#ifdef LBFGS
--- /dev/null
+ subroutine readpdb
+C Read the PDB file and convert the peptide geometry into virtual-chain
+C geometry.
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.NAMES'
+ include 'COMMON.CONTROL'
+ include 'COMMON.FRAG'
+ include 'COMMON.SETUP'
+ include 'COMMON.SBRIDGE'
+ character*3 seq,atom,res
+ character*80 card
+ double precision sccor(3,50)
+ double precision e1(3),e2(3),e3(3)
+ integer rescode,iterter(maxres),cou
+ logical fail
+ integer i,j,iii,ires,ires_old,ishift,ishift1,ibeg
+ double precision dcj
+ bfac=0.0d0
+ do i=1,maxres
+ iterter(i)=0
+ enddo
+ ibeg=1
+ ishift1=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
+ do
+ read (ipdbin,'(a80)',end=10) card
+ if (card(:5).eq.'HELIX') then
+ nhfrag=nhfrag+1
+ lsecondary=.true.
+ read(card(22:25),*) hfrag(1,nhfrag)
+ read(card(34:37),*) hfrag(2,nhfrag)
+ endif
+ if (card(:5).eq.'SHEET') then
+ nbfrag=nbfrag+1
+ lsecondary=.true.
+ read(card(24:26),*) bfrag(1,nbfrag)
+ read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc to be corrected !!!
+ bfrag(3,nbfrag)=bfrag(1,nbfrag)
+ bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
+ endif
+ if (card(:3).eq.'END') then
+ goto 10
+ else if (card(:3).eq.'TER') then
+C End current chain
+ ires_old=ires+2
+ itype(ires_old-1)=ntyp1
+ iterter(ires_old-1)=1
+ itype(ires_old)=ntyp1
+ iterter(ires_old)=1
+ ibeg=2
+ write (iout,*) "Chain ended",ires,ishift,ires_old
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+C Fish out the ATOM cards.
+c if (index(card(1:4),'ATOM').gt.0) then
+c read (card(14:16),'(a3)') atom
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+C Calculate the CM of the preceding residue.
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+ if (ires-ishift+ishift1.ne.ires_old) then
+ if (ibeg.eq.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires+nres)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+C Start new residue.
+c write (iout,'(a80)') card
+ if (ibeg.eq.1) then
+ ishift=ires-1
+ if (res.ne.'GLY' .and. res.ne. 'ACE') then
+ ishift=ishift-1
+ itype(1)=ntyp1
+ endif
+c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ibeg=0
+ else if (ibeg.eq.2) then
+c Start a new chain
+ ishift=-ires_old+ires-1
+c write (iout,*) "New chain started",ires,ishift
+ ibeg=0
+ endif
+ else
+ ires=ires-ishift
+c write (2,*) "ires",ires," ishift",ish
+ endif
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
+ if (res.eq.'ACE') then
+ itype(ires)=10
+ else
+ itype(ires)=rescode(ires,res,0)
+ endif
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+ read(card(61:66),*) bfac(ires)
+c if(me.eq.king.or..not.out1file)
+c & write (iout,'(2i3,2x,a,3f8.3)')
+c & ires,itype(ires),res,(c(j,ires),j=1,3)
+ iii=1
+ do j=1,3
+ sccor(j,iii)=c(j,ires)
+ enddo
+ else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N ' .and. atom.ne.'C ') then
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+ enddo
+ 10 if(me.eq.king.or..not.out1file)
+ & write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate dummy residue coordinates inside the "chain" of a multichain
+C system
+ nres=ires
+ do i=2,nres-1
+c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
+ if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
+ if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ print *,i,'tu dochodze'
+ call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif !fail
+ print *,i,'a tu?'
+ do j=1,3
+ c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i-2)-c(j,i-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i-1)+dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ else !itype(i+1).eq.ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i+3)-c(j,i+2))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i+1)-dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ endif !itype(i+1).eq.ntyp1
+ endif !itype.eq.ntyp1
+ enddo
+ write (iout,*) "After loop in readpbd"
+C Calculate the CM of the last side chain.
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ nsup=nres
+ nstart_sup=1
+ if (itype(nres).ne.10) then
+ nres=nres+1
+ itype(nres)=ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ endif
+ do i=2,nres-1
+ do j=1,3
+ c(j,i+nres)=dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ if (itype(1).eq.ntyp1) then
+ nsup=nsup-1
+ nstart_sup=2
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,4)-c(j,3))/2.0
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ endif
+C Calculate internal coordinates.
+ if(me.eq.king.or..not.out1file)then
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ endif
+ call flush(iout)
+c write(iout,*)"before int_from_cart nres",nres
+ call int_from_cart(.true.,.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
+c & vbld_inv(i+1)
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+c & vbld_inv(i+nres)
+ enddo
+ call sc_loc_geom(.false.)
+ call int_from_cart1(.false.)
+c call chainbuild
+C Copy the coordinates to reference coordinates
+ do i=1,nres
+ do j=1,3
+ cref(j,i)=c(j,i)
+ cref(j,i+nres)=c(j,i+nres)
+ enddo
+ enddo
+ 100 format (//' alpha-carbon coordinates ',
+ & ' centroid coordinates'/
+ 1 ' ', 6X,'X',11X,'Y',11X,'Z',
+ & 10X,'X',11X,'Y',11X,'Z')
+ 110 format (a,'(',i3,')',6f12.5)
+cc enddiag
+ do j=1,nbfrag
+ do i=1,4
+ bfrag(i,j)=bfrag(i,j)-ishift
+ enddo
+ enddo
+
+ do j=1,nhfrag
+ do i=1,2
+ hfrag(i,j)=hfrag(i,j)-ishift
+ enddo
+ enddo
+ return
+ end
+c---------------------------------------------------------------------------
+ subroutine readpdb_template(k)
+C Read the PDB file for read_constr_homology with read2sigma
+C and convert the peptide geometry into virtual-chain geometry.
+ implicit none
+ include 'DIMENSIONS'
+ include 'COMMON.LOCAL'
+ include 'COMMON.VAR'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.NAMES'
+ include 'COMMON.CONTROL'
+ include 'COMMON.FRAG'
+ include 'COMMON.SETUP'
+ integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+ & ishift_pdb,ires_ca
+ logical lprn /.false./,fail
+ double precision e1(3),e2(3),e3(3)
+ double precision dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
+ character*80 card
+ double precision sccor(3,20)
+ integer rescode,iterter(maxres)
+ do i=1,maxres
+ iterter(i)=0
+ enddo
+ ibeg=1
+ ishift1=0
+ ishift=0
+c write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
+ do
+ read (ipdbin,'(a80)',end=10) card
+ if (card(:3).eq.'END') then
+ goto 10
+ else if (card(:3).eq.'TER') then
+C End current chain
+ ires_old=ires+2
+ itype(ires_old-1)=ntyp1
+ iterter(ires_old-1)=1
+ itype(ires_old)=ntyp1
+ iterter(ires_old)=1
+ ibeg=2
+c write (iout,*) "Chain ended",ires,ishift,ires_old
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+C Fish out the ATOM cards.
+ if (index(card(1:4),'ATOM').gt.0) then
+ read (card(12:16),*) atom
+c write (iout,*) "! ",atom," !",ires
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+c write (iout,*) "ires",ires,ires-ishift+ishift1,
+c & " ires_old",ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
+C Calculate the CM of the preceding residue.
+ if (ibeg.eq.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
+ endif
+C Start new residue.
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+c write (iout,*) "BEG ires",ires
+ ishift=ires-1
+ if (res.ne.'GLY' .and. res.ne. 'ACE') then
+ ishift=ishift-1
+ itype(1)=ntyp1
+ endif
+ ires=ires-ishift+ishift1
+ ires_old=ires
+c write (iout,*) "ishift",ishift," ires",ires,
+c & " ires_old",ires_old
+c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ibeg=0
+ else if (ibeg.eq.2) then
+c Start a new chain
+ ishift=-ires_old+ires-1
+ ires=ires_old+1
+c write (iout,*) "New chain started",ires,ishift
+ ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
+ endif
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires)=10
+ else
+ itype(ires)=rescode(ires,res,0)
+ endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+c write (iout,*) "ires_old",ires_old," ires",ires
+c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c ishift1=ishift1+1
+c endif
+c write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)')
+ & ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
+ do j=1,3
+ sccor(j,iii)=c(j,ires)
+ enddo
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+c write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N' .and. atom.ne.'C' .and.
+ & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
+ & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c write (iout,*) "sidechain ",atom
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+ endif
+ enddo
+ 10 if(me.eq.king.or..not.out1file)
+ & write (iout,'(a,i5)') ' Nres: ',ires
+C Calculate dummy residue coordinates inside the "chain" of a multichain
+C system
+ nres=ires
+ do i=2,nres-1
+c write (iout,*) i,itype(i),itype(i+1)
+ if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
+ if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
+C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif !fail
+ do j=1,3
+ c(j,i)=c(j,i-1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i-2)-c(j,i-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i-1)+dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ else !itype(i+1).eq.ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i+3)-c(j,i+2))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i+1)-dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ endif !itype(i+1).eq.ntyp1
+ endif !itype.eq.ntyp1
+ enddo
+C Calculate the CM of the last side chain.
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ nsup=nres
+ nstart_sup=1
+ if (itype(nres).ne.10) then
+ nres=nres+1
+ itype(nres)=ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ endif
+ do i=2,nres-1
+ do j=1,3
+ c(j,i+nres)=dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ if (itype(1).eq.ntyp1) then
+ nsup=nsup-1
+ nstart_sup=2
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-1.9d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=(c(j,4)-c(j,3))/2.0
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ endif
+C Copy the coordinates to reference coordinates
+c do i=1,2*nres
+c do j=1,3
+c cref(j,i)=c(j,i)
+c enddo
+c enddo
+C Calculate internal coordinates.
+ if (out_template_coord) then
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+C Calculate internal coordinates.
+ call int_from_cart(.true.,.true.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
+c & vbld_inv(i+nres)
+ enddo
+ do i=1,nres
+ do j=1,3
+ cref(j,i)=c(j,i)
+ cref(j,i+nres)=c(j,i+nres)
+ enddo
+ enddo
+ do i=1,2*nres
+ do j=1,3
+ chomo(j,i,k)=c(j,i)
+ enddo
+ enddo
+
+ return
+ end
+
subroutine readpdb
C Read the PDB file and convert the peptide geometry into virtual-chain
C geometry.
- implicit none
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
- include 'COMMON.FRAG'
+ include 'COMMON.DISTFIT'
include 'COMMON.SETUP'
- include 'COMMON.SBRIDGE'
- character*3 seq,atom,res
+ include 'COMMON.FRAG'
+ integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+ & ishift_pdb
+ logical lprn /.false./,fail
+ double precision e1(3),e2(3),e3(3)
+ double precision dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
character*80 card
double precision sccor(3,50)
- double precision e1(3),e2(3),e3(3)
- integer rescode,iterter(maxres),cou
- logical fail
- integer i,j,iii,ires,ires_old,ishift,ibeg
- double precision dcj
- bfac=0.0d0
- do i=1,maxres
- iterter(i)=0
- enddo
+ integer rescode
+ efree_temp=0.0d0
ibeg=1
+ ishift1=0
+ ishift=0
+c write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
lsecondary=.false.
nhfrag=0
nbfrag=0
do
read (ipdbin,'(a80)',end=10) card
+c write (iout,'(a)') card
if (card(:5).eq.'HELIX') then
nhfrag=nhfrag+1
lsecondary=.true.
bfrag(4,nbfrag)=bfrag(2,nbfrag)
crc----------------------------------------
endif
- if (card(:3).eq.'END') then
- goto 10
- else if (card(:3).eq.'TER') then
-C End current chain
- ires_old=ires+2
- itype(ires_old-1)=ntyp1
- iterter(ires_old-1)=1
- itype(ires_old)=ntyp1
- iterter(ires_old)=1
- ibeg=2
- write (iout,*) "Chain ended",ires,ishift,ires_old
- if (unres_pdb) then
- do j=1,3
- dc(j,ires)=sccor(j,iii)
- enddo
- else
- call sccenter(ires,iii,sccor)
- endif
- endif
+ if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+c Read free energy
+ if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
- read (card(14:16),'(a3)') atom
- if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(12:16),*) atom
+c write (iout,*) "! ",atom," !",ires
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+c write (iout,*) "ires",ires,ires-ishift+ishift1,
+c & " ires_old",ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
+c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
if (ibeg.eq.0) then
+c write (iout,*) "Calculating sidechain center iii",iii
if (unres_pdb) then
do j=1,3
- dc(j,ires+nres)=sccor(j,iii)
+ dc(j,ires)=sccor(j,iii)
enddo
else
- call sccenter(ires,iii,sccor)
+ call sccenter(ires_old,iii,sccor)
endif
+ iii=0
endif
C Start new residue.
-c write (iout,'(a80)') card
- read (card(23:26),*) ires
- read (card(18:20),'(a3)') res
- if (ibeg.eq.1) then
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+c write (iout,*) "BEG ires",ires
ishift=ires-1
if (res.ne.'GLY' .and. res.ne. 'ACE') then
ishift=ishift-1
itype(1)=ntyp1
endif
-c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ires=ires-ishift+ishift1
+ ires_old=ires
+c write (iout,*) "ishift",ishift," ires",ires,
+c & " ires_old",ires_old
ibeg=0
- else if (ibeg.eq.2) then
-c Start a new chain
- ishift=-ires_old+ires-1
-c write (iout,*) "New chain started",ires,ishift
- ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
endif
- ires=ires-ishift
-c write (2,*) "ires",ires," ishift",ishift
- if (res.eq.'ACE') then
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
itype(ires)=10
else
itype(ires)=rescode(ires,res,0)
endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+c write (iout,*) "ires_old",ires_old," ires",ires
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c ishift1=ishift1+1
+ endif
+c write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
- read(card(61:66),*) bfac(ires)
-c if(me.eq.king.or..not.out1file)
-c & write (iout,'(2i3,2x,a,3f8.3)')
-c & ires,itype(ires),res,(c(j,ires),j=1,3)
- iii=1
+c write (iout,*) "backbone ",atom
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)')
+ & ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
do j=1,3
sccor(j,iii)=c(j,ires)
enddo
- else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
- & atom.ne.'N ' .and. atom.ne.'C ') then
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+c write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N' .and. atom.ne.'C' .and.
+ & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
+ & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c write (iout,*) "sidechain ",atom
iii=iii+1
read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
endif
endif
enddo
- 10 if(me.eq.king.or..not.out1file)
- & write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
- nres=ires
- do i=2,nres-1
-c write (iout,*) i,itype(i),itype(i+1),ntyp1,iterter(i)
- if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
- if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- print *,i,'tu dochodze'
- call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif !fail
- print *,i,'a tu?'
- do j=1,3
- c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i-2)-c(j,i-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i-1)+dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- else !itype(i+1).eq.ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
- call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i+3)-c(j,i+2))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i+1)-dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
- endif !itype.eq.ntyp1
- enddo
- write (iout,*) "After loop in readpbd"
+ 10 continue
+#ifdef DEBUG
+ write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+ if (ires.eq.0) return
C Calculate the CM of the last side chain.
+ if (iii.gt.0) then
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
enddo
- else
+ else
call sccenter(ires,iii,sccor)
endif
+ endif
+ nres=ires
nsup=nres
nstart_sup=1
if (itype(nres).ne.10) then
e2(3)=0.0d0
endif
do j=1,3
- c(j,nres)=c(j,nres-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
+ c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
enddo
else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
+ dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/dsqrt(2.0d0)
+ c(j,1)=c(j,2)-3.8d0*e2(j)
enddo
else
do j=1,3
- dcj=(c(j,4)-c(j,3))/2.0
+ dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
endif
endif
+C Copy the coordinates to reference coordinates
+c do i=1,2*nres
+c do j=1,3
+c cref(j,i)=c(j,i)
+c enddo
+c enddo
+C Calculate internal coordinates.
+ if (lprn) then
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ endif
C Calculate internal coordinates.
if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
do ires=1,nres
write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
& ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
& (c(j,nres+ires),j=1,3)
enddo
endif
- call flush(iout)
-c write(iout,*)"before int_from_cart nres",nres
call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
dc(j,i)=c(j,i+1)-c(j,i)
dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
enddo
-c write (iout,*) i,(dc(j,i),j=1,3),(dc_norm(j,i),j=1,3),
-c & vbld_inv(i+1)
enddo
do i=2,nres-1
do j=1,3
c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
c & vbld_inv(i+nres)
enddo
- call sc_loc_geom(.false.)
- call int_from_cart1(.false.)
c call chainbuild
C Copy the coordinates to reference coordinates
- do i=1,nres
+ do i=1,2*nres
do j=1,3
cref(j,i)=c(j,i)
- cref(j,i+nres)=c(j,i+nres)
enddo
enddo
- 100 format (//' alpha-carbon coordinates ',
- & ' centroid coordinates'/
- 1 ' ', 6X,'X',11X,'Y',11X,'Z',
- & 10X,'X',11X,'Y',11X,'Z')
- 110 format (a,'(',i3,')',6f12.5)
-cc enddiag
+
+
do j=1,nbfrag
do i=1,4
bfrag(i,j)=bfrag(i,j)-ishift
hfrag(i,j)=hfrag(i,j)-ishift
enddo
enddo
+ ishift_pdb=ishift
return
end
-c---------------------------------------------------------------------------
- subroutine int_from_cart(lside,lprn)
- implicit none
- include 'DIMENSIONS'
-#ifdef MPI
- include "mpif.h"
-#endif
- include 'COMMON.LOCAL'
- include 'COMMON.VAR'
- include 'COMMON.CHAIN'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.GEO'
- include 'COMMON.NAMES'
- include 'COMMON.CONTROL'
- include 'COMMON.SETUP'
- double precision dist,alpha,beta
- character*3 seq,atom,res
- character*80 card
- double precision sccor(3,50)
- integer rescode
- logical lside,lprn
- integer i,j,iti
- double precision di,cosfac2,sinfac2,cosfac,sinfac
-#ifdef MPI
- if(me.eq.king.or..not.out1file)then
-#endif
- if (lprn) then
- write (iout,'(/a)')
- & 'Internal coordinates calculated from crystal structure.'
- if (lside) then
- write (iout,'(8a)') ' Res ',' dvb',' Theta',
- & ' Phi',' Dsc_id',' Dsc',' Alpha',
- & ' Omega'
- else
- write (iout,'(4a)') ' Res ',' dvb',' Theta',
- & ' Phi'
- endif
- endif
-#ifdef MPI
- endif
-#endif
- do i=1,nres-1
- iti=itype(i)
- if (iti.ne.ntyp1 .and. itype(i+1).ne.ntyp1 .and.
- & (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0)) then
- write (iout,'(a,i4)') 'Bad Cartesians for residue',i
-ctest stop
- endif
- vbld(i+1)=dist(i,i+1)
- vbld_inv(i+1)=1.0d0/vbld(i+1)
-c write (iout,*) "i",i+1," vbld",vbld(i+1)," vbld_inv",
-c & vbld_inv(i+1)
- if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
- if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
- enddo
-c if (unres_pdb) then
-c if (itype(1).eq.21) then
-c theta(3)=90.0d0*deg2rad
-c phi(4)=180.0d0*deg2rad
-c vbld(2)=3.8d0
-c vbld_inv(2)=1.0d0/vbld(2)
-c endif
-c if (itype(nres).eq.21) then
-c theta(nres)=90.0d0*deg2rad
-c phi(nres)=180.0d0*deg2rad
-c vbld(nres)=3.8d0
-c vbld_inv(nres)=1.0d0/vbld(2)
-c endif
-c endif
-c print *,"A TU2"
- if (lside) then
- do i=2,nres-1
- do j=1,3
- c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
- & +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
- enddo
- iti=itype(i)
- di=dist(i,nres+i)
- vbld(i+nres)=di
- if (itype(i).ne.10) then
- vbld_inv(i+nres)=1.0d0/di
- else
- vbld_inv(i+nres)=0.0d0
- endif
- if (iti.ne.10) then
- alph(i)=alpha(nres+i,i,maxres2)
- omeg(i)=beta(nres+i,i,maxres2,i+1)
- endif
- if(me.eq.king.or..not.out1file)then
- if (lprn)
- & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
- & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
- & rad2deg*alph(i),rad2deg*omeg(i)
- endif
- enddo
- if (lprn) then
- i=nres
- iti=itype(nres)
- write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
- & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
- & rad2deg*alph(i),rad2deg*omeg(i)
- endif
- else if (lprn) then
- do i=2,nres
- iti=itype(i)
- if(me.eq.king.or..not.out1file)
- & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
- & rad2deg*theta(i),rad2deg*phi(i)
- enddo
- endif
- return
- end
-c-------------------------------------------------------------------------------
- subroutine sc_loc_geom(lprn)
- implicit none
- include 'DIMENSIONS'
-#ifdef MPI
- include "mpif.h"
-#endif
- include 'COMMON.LOCAL'
- include 'COMMON.VAR'
- include 'COMMON.CHAIN'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.GEO'
- include 'COMMON.NAMES'
- include 'COMMON.CONTROL'
- include 'COMMON.SETUP'
- double precision x_prime(3),y_prime(3),z_prime(3)
- logical lprn
- integer i,j,it
- double precision xx,yy,zz,cosfac,cosfac2,sinfac,sinfac2
- do i=1,nres-1
- do j=1,3
- dc_norm(j,i)=vbld_inv(i+1)*(c(j,i+1)-c(j,i))
- enddo
-c write (iout,*) "i",i," dc",(dc_norm(j,i),j=1,3),
-c & " vbld",vbld_inv(i+1)
- enddo
- do i=2,nres-1
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
- do j=1,3
- dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
- enddo
-c write (iout,*) "i",i," dc",(dc_norm(j,i+nres),j=1,3),
-c & " vbld",vbld_inv(i+nres)
- else
- do j=1,3
- dc_norm(j,i+nres)=0.0d0
- enddo
- endif
- enddo
- do i=2,nres-1
- costtab(i+1) =dcos(theta(i+1))
- sinttab(i+1) =dsqrt(1-costtab(i+1)*costtab(i+1))
- cost2tab(i+1)=dsqrt(0.5d0*(1.0d0+costtab(i+1)))
- sint2tab(i+1)=dsqrt(0.5d0*(1.0d0-costtab(i+1)))
- cosfac2=0.5d0/(1.0d0+costtab(i+1))
- cosfac=dsqrt(cosfac2)
- sinfac2=0.5d0/(1.0d0-costtab(i+1))
- sinfac=dsqrt(sinfac2)
- it=itype(i)
-c write (iout,*) "i",i," costab",costtab(i+1),
-c & " sintab",sinttab(i+1)
-c write (iout,*) "dc_norm_b",(dc_norm(j,i-1),j=1,3)
-c write (iout,*) "dc_norm_s",(dc_norm(j,i+nres),j=1,3)
- if (it.ne.10 .and. itype(i).ne.ntyp1) then
-c
-C Compute the axes of tghe local cartesian coordinates system; store in
-c x_prime, y_prime and z_prime
-c
- do j=1,3
- x_prime(j) = 0.00
- y_prime(j) = 0.00
- z_prime(j) = 0.00
- enddo
- do j = 1,3
- x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
- y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
- enddo
-c write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c write (iout,*) "y_prime",(y_prime(j),j=1,3)
- call vecpr(x_prime,y_prime,z_prime)
-c write (iout,*) "z_prime",(z_prime(j),j=1,3)
-c
-C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
-C to local coordinate system. Store in xx, yy, zz.
-c
- xx=0.0d0
- yy=0.0d0
- zz=0.0d0
- do j = 1,3
- xx = xx + x_prime(j)*dc_norm(j,i+nres)
- yy = yy + y_prime(j)*dc_norm(j,i+nres)
- zz = zz + z_prime(j)*dc_norm(j,i+nres)
- enddo
-
- xxref(i)=xx
- yyref(i)=yy
- zzref(i)=zz
- else
- xxref(i)=0.0d0
- yyref(i)=0.0d0
- zzref(i)=0.0d0
- endif
- enddo
- if (lprn) then
-#ifdef MPI
- if (me.eq.king.or..not.out1file) then
-#endif
- write (iout,*) "xxref,yyref,zzref"
- do i=2,nres
- write (iout,'(a3,i4,3f10.5)')
- & restyp(itype(i)),i,xxref(i),yyref(i),zzref(i)
- enddo
-#ifdef MPI
- endif
-#endif
- endif
- return
- end
-c---------------------------------------------------------------------------
- subroutine sccenter(ires,nscat,sccor)
- implicit none
- include 'DIMENSIONS'
- include 'COMMON.CHAIN'
- integer i,j,ires,nscat
- double precision sccor(3,50)
- double precision sccmj
- do j=1,3
- sccmj=0.0D0
- do i=1,nscat
- sccmj=sccmj+sccor(j,i)
- enddo
- dc(j,ires)=sccmj/nscat
- enddo
- return
- end
-c---------------------------------------------------------------------------
- subroutine bond_regular
- implicit none
- include 'DIMENSIONS'
- include 'COMMON.VAR'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.CHAIN'
- integer i,i1,i2
- do i=1,nres-1
- vbld(i+1)=vbl
- vbld_inv(i+1)=vblinv
- vbld(i+1+nres)=dsc(iabs(itype(i+1)))
- vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
-c print *,vbld(i+1),vbld(i+1+nres)
- enddo
-c Adam 2/26/20 Alter virtual bonds for non-blocking end groups of each chain
- do i=1,nchain
- i1=chain_border(1,i)
- i2=chain_border(2,i)
- if (i1.gt.1) then
- vbld(i1)=vbld(i1)/2
- vbld_inv(i1)=vbld_inv(i1)*2
- endif
- if (i2.lt.nres) then
- vbld(i2+1)=vbld(i2+1)/2
- vbld_inv(i2+1)=vbld_inv(i2+1)*2
- endif
- enddo
- return
- end
-c---------------------------------------------------------------------------
+c-----------------------------------------------------------------------
subroutine readpdb_template(k)
-C Read the PDB file for read_constr_homology with read2sigma
+C Read the PDB file with gaps for read_constr_homology with read2sigma
C and convert the peptide geometry into virtual-chain geometry.
- implicit none
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.GEO'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
- include 'COMMON.FRAG'
+ include 'COMMON.DISTFIT'
include 'COMMON.SETUP'
- integer i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
- & ishift_pdb,ires_ca
+ include 'COMMON.MD'
+ integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity,
+ & ishift_pdb
logical lprn /.false./,fail
double precision e1(3),e2(3),e3(3)
double precision dcj,efree_temp
character*3 seq,res
character*5 atom
character*80 card
- double precision sccor(3,20)
- integer rescode,iterter(maxres)
- do i=1,maxres
- iterter(i)=0
- enddo
+ double precision sccor(3,50)
+ integer rescode
+ efree_temp=0.0d0
ibeg=1
ishift1=0
ishift=0
lsecondary=.false.
nhfrag=0
nbfrag=0
- do
+ do
read (ipdbin,'(a80)',end=10) card
- if (card(:3).eq.'END') then
- goto 10
- else if (card(:3).eq.'TER') then
-C End current chain
- ires_old=ires+2
- itype(ires_old-1)=ntyp1
- iterter(ires_old-1)=1
- itype(ires_old)=ntyp1
- iterter(ires_old)=1
- ibeg=2
-c write (iout,*) "Chain ended",ires,ishift,ires_old
- if (unres_pdb) then
- do j=1,3
- dc(j,ires)=sccor(j,iii)
- enddo
- else
- call sccenter(ires,iii,sccor)
- endif
- endif
+c write (iout,'(a)') card
+ if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
read (card(12:16),*) atom
c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
+c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
if (ibeg.eq.0) then
+c write (iout,*) "Calculating sidechain center iii",iii
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
ires_old=ires
c write (iout,*) "ishift",ishift," ires",ires,
c & " ires_old",ires_old
-c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
- ibeg=0
- else if (ibeg.eq.2) then
-c Start a new chain
- ishift=-ires_old+ires-1
- ires=ires_old+1
-c write (iout,*) "New chain started",ires,ishift
ibeg=0
else
ishift=ishift-(ires-ishift+ishift1-ires_old-1)
ires=ires-ishift+ishift1
endif
c write (iout,*) "ires_old",ires_old," ires",ires
-c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
c ishift1=ishift1+1
-c endif
+ endif
c write (2,*) "ires",ires," res ",res," ity",ity
if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
& res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+c write (iout,*) "backbone ",atom
#ifdef DEBUG
write (iout,'(2i3,2x,a,3f8.3)')
& ires,itype(ires),res,(c(j,ires),j=1,3)
endif
endif
enddo
- 10 if(me.eq.king.or..not.out1file)
- & write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
- nres=ires
- do i=2,nres-1
-c write (iout,*) i,itype(i),itype(i+1)
- if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
- if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif !fail
- do j=1,3
- c(j,i)=c(j,i-1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i-2)-c(j,i-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i-1)+dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- else !itype(i+1).eq.ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
- call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i+3)-c(j,i+2))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i+1)-dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
- endif !itype.eq.ntyp1
- enddo
+ 10 continue
+#ifdef DEBUG
+ write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+ if (ires.eq.0) return
C Calculate the CM of the last side chain.
+ if (iii.gt.0) then
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
else
call sccenter(ires,iii,sccor)
endif
+ endif
+ nres=ires
nsup=nres
nstart_sup=1
if (itype(nres).ne.10) then
nres=nres+1
itype(nres)=ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
- enddo
- else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
+ dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
endif
- endif
do i=2,nres-1
do j=1,3
c(j,i+nres)=dc(j,i)
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-1.9d0*e2(j)
+ c(j,1)=c(j,2)-3.8d0*e2(j)
enddo
else
do j=1,3
- dcj=(c(j,4)-c(j,3))/2.0
+ dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
endif
endif
-C Copy the coordinates to reference coordinates
-c do i=1,2*nres
-c do j=1,3
-c cref(j,i)=c(j,i)
-c enddo
-c enddo
C Calculate internal coordinates.
- if (out_template_coord) then
+ if (lprn) then
write (iout,'(/a)')
& "Cartesian coordinates of the reference structure"
write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
enddo
endif
C Calculate internal coordinates.
- call int_from_cart(.true.,.true.)
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ endif
+ call int_from_cart(.true.,.false.)
call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
c & vbld_inv(i+nres)
enddo
- do i=1,nres
- do j=1,3
- cref(j,i)=c(j,i)
- cref(j,i+nres)=c(j,i+nres)
- enddo
- enddo
+c call chainbuild
+C Copy the coordinates to reference coordinates
do i=1,2*nres
do j=1,3
+ cref(j,i)=c(j,i)
chomo(j,i,k)=c(j,i)
enddo
enddo
+
+ ishift_pdb=ishift
return
end
-
dccart=(index(controlcard,'CART').gt.0)
overlapsc=(index(controlcard,'OVERLAP').gt.0)
overlapsc=.not.overlapsc
- searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
- searchsc=.not.searchsc
+ searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
+ & index(controlcard,'NOSEARCHSC').eq.0)
sideadd=(index(controlcard,'SIDEADD').gt.0)
energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
if (iz_sc.lt.2) then
do ichain=1,nchain
indchain=tabpermchain(ichain,iperm)
+#define DEBUG
#ifdef DEBUG
write (iout,*) "ichain",ichain," indchain",indchain
write (iout,*) "chain_border",chain_border(1,ichain),
enddo
rmsside=dsqrt(rmsside/nnnn)
#ifdef DEBUG
- write (iout,*) iii,iref," nnnn",nnnn," rmsside",rmsside
+ write (iout,*) "nnnn",nnnn," rmsside",rmsside
#endif
+#undef DEBUG
rmscalc_side=rmsside
return
end
mask_side=1
crc n_fun=n_fun+1
ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime
+ call chainbuild_extconf
return
end
integer IERROR
integer i,j,k,l,ind,ind1,m,ii,iti,it,nzero,innt,inct
integer ichain,nind
- logical lprn /.true./
+ logical lprn /.false./
double precision dtdi,gamvec(MAXRES2)
common /syfek/ gamvec
#ifndef FIVEDIAG
endif
if (overlapsc) then
-c print *, 'Calling OVERLAP_SC'
+ write (iout,*) 'Calling OVERLAP_SC'
call overlap_sc(fail)
-c print *,"After overlap_sc"
+ write (iout,*) "After overlap_sc"
endif
if (searchsc) then
BIN = ~/bin
FC = ftn
-OPT = -mcmodel=medium -shared-intel -O3 -dynamic
+#OPT = -mcmodel=medium -shared-intel -O3 -dynamic
#OPT = -O3 -intel-static -mcmodel=medium
#OPT = -O3 -ip -w
-#OPT = -g -CB -mcmodel=medium -shared-intel -dynamic
+OPT = -g -CB -mcmodel=medium -shared-intel -dynamic
FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include
LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a
subroutine readpdb
C Read the PDB file and convert the peptide geometry into virtual-chain
C geometry.
- implicit none
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'DIMENSIONS.ZSCOPT'
- include 'COMMON.CONTROL'
+ include 'DIMENSIONS.FREE'
+ include 'COMMON.FRAG'
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
- include 'COMMON.SBRIDGE'
- character*3 seq,atom,res
+ include 'COMMON.CONTROL'
+ integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
+ logical lprn /.false./,fail
+ double precision e1(3),e2(3),e3(3)
+ double precision dcj,efree_temp
+ character*3 seq,res
+ character*5 atom
character*80 card
- double precision sccor(3,50)
- integer i,j,iii,ibeg,ishift,ishift1,ity,ires,ires_old
- double precision dcj
- integer rescode,kkk,lll,icha,cou,kupa,iprzes
+ double precision sccor(3,20)
+ integer rescode
+ efree_temp=0.0d0
ibeg=1
ishift1=0
+ ishift=0
+c write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
do
read (ipdbin,'(a80)',end=10) card
- if (card(:3).eq.'END') then
- goto 10
- else if (card(:3).eq.'TER') then
-C End current chain
-c ires_old=ires+1
- ires_old=ires+2
- itype(ires_old-1)=ntyp1
- itype(ires_old)=ntyp1
- ibeg=2
-c write (iout,*) "Chain ended",ires,ishift,ires_old
- call sccenter(ires,iii,sccor)
+c write (iout,'(a)') card
+ if (card(:5).eq.'HELIX') then
+ nhfrag=nhfrag+1
+ lsecondary=.true.
+ read(card(22:25),*) hfrag(1,nhfrag)
+ read(card(34:37),*) hfrag(2,nhfrag)
+ endif
+ if (card(:5).eq.'SHEET') then
+ nbfrag=nbfrag+1
+ lsecondary=.true.
+ read(card(24:26),*) bfrag(1,nbfrag)
+ read(card(35:37),*) bfrag(2,nbfrag)
+crc----------------------------------------
+crc to be corrected !!!
+ bfrag(3,nbfrag)=bfrag(1,nbfrag)
+ bfrag(4,nbfrag)=bfrag(2,nbfrag)
+crc----------------------------------------
endif
+ if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
+c Read free energy
+ if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
- read (card(14:16),'(a3)') atom
- if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(12:16),*) atom
+c write (iout,*) "! ",atom," !",ires
+c if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+c write (iout,*) "ires",ires,ires-ishift+ishift1,
+c & " ires_old",ires_old
+c write (iout,*) "ishift",ishift," ishift1",ishift1
+c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
+c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
if (ibeg.eq.0) then
- call sccenter(ires,iii,sccor)
+c write (iout,*) "Calculating sidechain center iii",iii
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
endif
C Start new residue.
-c write (iout,'(a80)') card
- read (card(23:26),*) ires
- read (card(18:20),'(a3)') res
- if (ibeg.eq.1) then
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+c write (iout,*) "BEG ires",ires
ishift=ires-1
if (res.ne.'GLY' .and. res.ne. 'ACE') then
ishift=ishift-1
itype(1)=ntyp1
endif
-c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
+ ires=ires-ishift+ishift1
+ ires_old=ires
+c write (iout,*) "ishift",ishift," ires",ires,
+c & " ires_old",ires_old
ibeg=0
- else if (ibeg.eq.2) then
-c Start a new chain
- ishift=-ires_old+ires-1
-c write (iout,*) "New chain started",ires,ishift
- ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
endif
- ires=ires-ishift
-c write (2,*) "ires",ires," ishift",ishift
- if (res.eq.'ACE') then
- ity=10
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires)=10
else
itype(ires)=rescode(ires,res,0)
endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+c write (iout,*) "ires_old",ires_old," ires",ires
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+c ishift1=ishift1+1
+ endif
+c write (2,*) "ires",ires," res ",res," ity",ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
+ & res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
- read(card(61:66),*) bfac(ires)
-c write (iout,'(2i3,2x,a,3f8.3,5x,f8.3)')
-c & ires,itype(ires),res,(c(j,ires),j=1,3),bfac(ires)
- iii=1
+c write (iout,*) "backbone ",atom
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)')
+ & ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
do j=1,3
sccor(j,iii)=c(j,ires)
enddo
- else if (atom.ne.'O '.and.atom(1:1).ne.'H' .and.
- & atom(1:1).ne.'Q' .and. atom(1:2).ne.'1H' .and.
- & atom(1:2).ne.'2H' .and. atom(1:2).ne.'3H' .and.
- & atom.ne.'N ' .and. atom.ne.'C ' .and.
- & atom.ne.'OXT' ) then
+ if (ishift.ne.0) then
+ ires_ca=ires+ishift-ishift1
+ else
+ ires_ca=ires
+ endif
+c write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.
+ & atom.ne.'N' .and. atom.ne.'C' .and.
+ & atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.
+ & atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+c write (iout,*) "sidechain ",atom
iii=iii+1
-c write (iout,*) res,ires,iii,atom
read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
-c write (iout,'(3f8.3)') (sccor(j,iii),j=1,3)
endif
endif
enddo
- 10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
- nres=ires
- do i=2,nres-1
-c write (iout,*) i,itype(i)
-
- if (itype(i).eq.ntyp1) then
- if (itype(i+1).eq.ntyp1) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
-C if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
-C call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
-C if (fail) then
-C e2(1)=0.0d0
-C e2(2)=1.0d0
-C e2(3)=0.0d0
-C endif !fail
-C do j=1,3
-C c(j,i)=c(j,i-1)-1.9d0*e2(j)
-C enddo
-C else !unres_pdb
- do j=1,3
- dcj=(c(j,i-2)-c(j,i-3))/2.0
- c(j,i)=c(j,i-1)+dcj
- c(j,nres+i)=c(j,i)
- enddo
-C endif !unres_pdb
- else !itype(i+1).eq.ntyp1
-C if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
-C call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
-C if (fail) then
-C e2(1)=0.0d0
-C e2(2)=1.0d0
-C e2(3)=0.0d0
-C endif
-C do j=1,3
-C c(j,i)=c(j,i+1)-1.9d0*e2(j)
-C enddo
-C else !unres_pdb
- do j=1,3
- dcj=(c(j,i+3)-c(j,i+2))/2.0
- c(j,i)=c(j,i+1)-dcj
- c(j,nres+i)=c(j,i)
- enddo
-C endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
- endif !itype.eq.ntyp1
- enddo
+ 10 continue
+#ifdef DEBUG
+ write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+ if (ires.eq.0) return
C Calculate the CM of the last side chain.
- call sccenter(ires,iii,sccor)
+ if (iii.gt.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+ nres=ires
nsup=nres
nstart_sup=1
if (itype(nres).ne.10) then
nres=nres+1
itype(nres)=ntyp1
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+ enddo
+ else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
+ dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
+ endif
endif
do i=2,nres-1
do j=1,3
if (itype(1).eq.ntyp1) then
nsup=nsup-1
nstart_sup=2
+ if (unres_pdb) then
+C 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-3.8d0*e2(j)
+ enddo
+ else
do j=1,3
- dcj=(c(j,4)-c(j,3))/2.0
+ dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
+ endif
endif
+C Copy the coordinates to reference coordinates
+c do i=1,2*nres
+c do j=1,3
+c cref(j,i)=c(j,i)
+c enddo
+c enddo
C Calculate internal coordinates.
- write (iout,100)
+ if (lprn) then
+ write (iout,'(/a)')
+ & "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
+ & "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)')
+ & restyp(itype(ires)),ires,(c(j,ires),j=1,3),
+ & (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+C Calculate internal coordinates.
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
& ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
& (c(j,nres+ires),j=1,3)
- enddo
+ enddo
call int_from_cart(.true.,.false.)
- call flush(iout)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
do i=1,nres-1
do j=1,3
dc(j,i)=c(j,i+1)-c(j,i)
enddo
c call chainbuild
C Copy the coordinates to reference coordinates
- do i=1,nres
+ do i=1,2*nres
do j=1,3
cref(j,i)=c(j,i)
- cref(j,i+nres)=c(j,i+nres)
enddo
enddo
- 100 format ('Residue alpha-carbon coordinates ',
- & ' centroid coordinates'/
- 1 ' ', 6X,'X',7X,'Y',7X,'Z',
- & 12X,'X',7X,'Y',7X,'Z')
- 110 format (a,'(',i3,')',6f12.5)
+
+ do j=1,nbfrag
+ do i=1,4
+ bfrag(i,j)=bfrag(i,j)-ishift
+ enddo
+ enddo
+
+ do j=1,nhfrag
+ do i=1,2
+ hfrag(i,j)=hfrag(i,j)-ishift
+ enddo
+ enddo
ishift_pdb=ishift
return
end
c---------------------------------------------------------------------------
subroutine int_from_cart(lside,lprn)
- implicit none
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'DIMENSIONS.ZSCOPT'
include 'COMMON.LOCAL'
include 'COMMON.IOUNITS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
- character*3 seq,atom,res
+ include 'COMMON.CONTROL'
+ character*3 seq,res
+c character*5 atom
character*80 card
- double precision sccor(3,50)
+ dimension sccor(3,20)
integer rescode
- double precision dist,alpha,beta,di
- integer i,j,iti
logical lside,lprn
- if (lprn) then
+ if (lprn) then
write (iout,'(/a)')
& 'Internal coordinates calculated from crystal structure.'
if (lside) then
write (iout,'(8a)') ' Res ',' dvb',' Theta',
- & ' Phi',' Dsc_id',' Dsc',' Alpha',
- & ' Omega'
+ & ' Gamma',' Dsc_id',' Dsc',' Alpha',
+ & ' Beta '
else
write (iout,'(4a)') ' Res ',' dvb',' Theta',
- & ' Phi'
+ & ' Gamma'
endif
- endif
- do i=2,nres
+ endif
+ do i=1,nres-1
iti=itype(i)
-c write (iout,*) i,i-1,(c(j,i),j=1,3),(c(j,i-1),j=1,3),dist(i,i-1)
- if (itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1 .and.
- & (dist(i,i-1).lt.1.0D0 .or. dist(i,i-1).gt.6.0D0)) then
+ if (dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0) then
write (iout,'(a,i4)') 'Bad Cartesians for residue',i
- stop
+ctest stop
endif
- vbld(i)=dist(i-1,i)
- vbld_inv(i)=1.0d0/vbld(i)
- theta(i+1)=alpha(i-1,i,i+1)
+ vbld(i+1)=dist(i,i+1)
+ vbld_inv(i+1)=1.0d0/vbld(i+1)
+ if (i.gt.1) theta(i+1)=alpha(i-1,i,i+1)
if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
enddo
-c if (itype(1).eq.ntyp1) then
-c do j=1,3
-c c(j,1)=c(j,2)+(c(j,3)-c(j,4))
-c enddo
-c endif
-c if (itype(nres).eq.ntyp1) then
-c do j=1,3
-c c(j,nres)=c(j,nres-1)+(c(j,nres-2)-c(j,nres-3))
-c enddo
+c if (unres_pdb) then
+c if (itype(1).eq.ntyp1) then
+c theta(3)=90.0d0*deg2rad
+c phi(4)=180.0d0*deg2rad
+c vbld(2)=3.8d0
+c vbld_inv(2)=1.0d0/vbld(2)
+c endif
+c if (itype(nres).eq.ntyp1) then
+c theta(nres)=90.0d0*deg2rad
+c phi(nres)=180.0d0*deg2rad
+c vbld(nres)=3.8d0
+c vbld_inv(nres)=1.0d0/vbld(2)
+c endif
c endif
if (lside) then
do i=2,nres-1
do j=1,3
- c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
+ c(j,maxres2)=0.5D0*(2*c(j,i)+(c(j,i-1)-c(j,i))*vbld_inv(i)
+ & +(c(j,i+1)-c(j,i))*vbld_inv(i+1))
enddo
iti=itype(i)
di=dist(i,nres+i)
- vbld(i+nres)=di
+C 10/03/12 Adam: Correction for zero SC-SC bond length
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1. and. di.eq.0.0d0)
+ & di=dsc(itype(i))
+ vbld(i+nres)=di
if (itype(i).ne.10) then
vbld_inv(i+nres)=1.0d0/di
else
alph(i)=alpha(nres+i,i,maxres2)
omeg(i)=beta(nres+i,i,maxres2,i+1)
endif
- if (iti.ne.10) then
- alph(i)=alpha(nres+i,i,maxres2)
- omeg(i)=beta(nres+i,i,maxres2,i+1)
- endif
- if (lprn)
- & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
- & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),di,
- & rad2deg*alph(i),rad2deg*omeg(i)
+ if (lprn)
+ & write (iout,'(a3,i4,7f10.3)') restyp(iti),i,vbld(i),
+ & rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),
+ & rad2deg*alph(i),rad2deg*omeg(i)
enddo
else if (lprn) then
do i=2,nres
iti=itype(i)
write (iout,'(a3,i4,7f10.3)') restyp(iti),i,dist(i,i-1),
- & rad2deg*theta(i),rad2deg*phi(i)
+ & rad2deg*theta(i),rad2deg*phi(i)
enddo
endif
return
end
-c---------------------------------------------------------------------------
- subroutine sccenter(ires,nscat,sccor)
- implicit none
- include 'DIMENSIONS'
- include 'COMMON.CHAIN'
- integer ires,nscat,i,j
- double precision sccor(3,50),sccmj
- do j=1,3
- sccmj=0.0D0
- do i=1,nscat
- sccmj=sccmj+sccor(j,i)
- enddo
- dc(j,ires)=sccmj/nscat
- enddo
- return
- end
-c---------------------------------------------------------------------------
+c-------------------------------------------------------------------------------
subroutine sc_loc_geom(lprn)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
- include 'COMMON.SETUP'
double precision x_prime(3),y_prime(3),z_prime(3)
logical lprn
do i=1,nres-1
enddo
enddo
do i=2,nres-1
- if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ if (itype(i).ne.10) then
do j=1,3
dc_norm(j,i+nres)=vbld_inv(i+nres)*(c(j,i+nres)-c(j,i))
enddo
sinfac2=0.5d0/(1.0d0-costtab(i+1))
sinfac=dsqrt(sinfac2)
it=itype(i)
- if (it.ne.10 .and. itype(i).ne.ntyp1) then
+ if (it.ne.10) then
c
C Compute the axes of tghe local cartesian coordinates system; store in
c x_prime, y_prime and z_prime
x_prime(j) = (dc_norm(j,i) - dc_norm(j,i-1))*cosfac
y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
enddo
-c write (iout,*) "x_prime",(x_prime(j),j=1,3)
-c write (iout,*) "y_prime",(y_prime(j),j=1,3)
call vecpr(x_prime,y_prime,z_prime)
-c write (iout,*) "z_prime",(z_prime(j),j=1,3)
c
C Transform the unit vector of the ith side-chain centroid, dC_norm(*,i),
C to local coordinate system. Store in xx, yy, zz.
endif
enddo
if (lprn) then
- write (iout,*) "xxref,yyref,zzref"
do i=2,nres
iti=itype(i)
- write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),yyref(i),
- & zzref(i)
+ write (iout,'(a3,i4,3f10.5)') restyp(iti),i,xxref(i),
+ & yyref(i),zzref(i)
enddo
endif
return
end
c---------------------------------------------------------------------------
+ subroutine sccenter(ires,nscat,sccor)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CHAIN'
+ dimension sccor(3,20)
+ do j=1,3
+ sccmj=0.0D0
+ do i=1,nscat
+ sccmj=sccmj+sccor(j,i)
+ enddo
+ dc(j,ires)=sccmj/nscat
+ enddo
+ return
+ end
+c---------------------------------------------------------------------------
subroutine bond_regular
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
do i=1,nres-1
vbld(i+1)=vbl
vbld_inv(i+1)=1.0d0/vbld(i+1)
- vbld(i+1+nres)=dsc(iabs(itype(i+1)))
- vbld_inv(i+1+nres)=dsc_inv(iabs(itype(i+1)))
+ vbld(i+1+nres)=dsc(itype(i+1))
+ vbld_inv(i+1+nres)=dsc_inv(itype(i+1))
c print *,vbld(i+1),vbld(i+1+nres)
enddo
return
end
-c---------------------------------------------------------------------------
subroutine readpdb_template(k)
-C Read the PDB file for read_constr_homology with read2sigma
+C Read the PDB file with gaps for read_constr_homology with read2sigma
C and convert the peptide geometry into virtual-chain geometry.
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
include 'COMMON.LOCAL'
include 'COMMON.VAR'
include 'COMMON.CHAIN'
include 'COMMON.GEO'
include 'COMMON.NAMES'
include 'COMMON.CONTROL'
- include 'COMMON.SETUP'
integer i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity
logical lprn /.false./,fail
double precision e1(3),e2(3),e3(3)
character*5 atom
character*80 card
double precision sccor(3,20)
- integer rescode,iterter(maxres)
- do i=1,maxres
- iterter(i)=0
- enddo
+ integer rescode
+ efree_temp=0.0d0
ibeg=1
ishift1=0
ishift=0
lsecondary=.false.
nhfrag=0
nbfrag=0
- do
+ do
read (ipdbin,'(a80)',end=10) card
- if (card(:3).eq.'END') then
- goto 10
- else if (card(:3).eq.'TER') then
-C End current chain
- ires_old=ires+2
- itype(ires_old-1)=ntyp1
- iterter(ires_old-1)=1
- itype(ires_old)=ntyp1
- iterter(ires_old)=1
- ibeg=2
-c write (iout,*) "Chain ended",ires,ishift,ires_old
- if (unres_pdb) then
- do j=1,3
- dc(j,ires)=sccor(j,iii)
- enddo
- else
- call sccenter(ires,iii,sccor)
- endif
- endif
+c write (iout,'(a)') card
+ if (card(:3).eq.'END' .or. card(:3).eq.'TER') goto 10
C Fish out the ATOM cards.
if (index(card(1:4),'ATOM').gt.0) then
read (card(12:16),*) atom
c write (iout,*) "IRES",ires-ishift+ishift1,ires_old
if (ires-ishift+ishift1.ne.ires_old) then
C Calculate the CM of the preceding residue.
+c if (ibeg.eq.0) call sccenter(ires,iii,sccor)
if (ibeg.eq.0) then
+c write (iout,*) "Calculating sidechain center iii",iii
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
ires_old=ires
c write (iout,*) "ishift",ishift," ires",ires,
c & " ires_old",ires_old
-c write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
- ibeg=0
- else if (ibeg.eq.2) then
-c Start a new chain
- ishift=-ires_old+ires-1
- ires=ires_old+1
-c write (iout,*) "New chain started",ires,ishift
ibeg=0
else
ishift=ishift-(ires-ishift+ishift1-ires_old-1)
ires=ires-ishift+ishift1
endif
c write (iout,*) "ires_old",ires_old," ires",ires
-c if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
c ishift1=ishift1+1
-c endif
+ endif
c write (2,*) "ires",ires," res ",res," ity",ity
if (atom.eq.'CA' .or. atom.eq.'CH3' .or.
& res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-c write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
+c write (iout,*) "backbone ",atom
#ifdef DEBUG
write (iout,'(2i3,2x,a,3f8.3)')
& ires,itype(ires),res,(c(j,ires),j=1,3)
endif
endif
enddo
- 10 write (iout,'(a,i5)') ' Nres: ',ires
-C Calculate dummy residue coordinates inside the "chain" of a multichain
-C system
- nres=ires
- do i=2,nres-1
-c write (iout,*) i,itype(i),itype(i+1)
- if (itype(i).eq.ntyp1.and.iterter(i).eq.1) then
- if (itype(i+1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
-C 16/01/2014 by Adasko: Adding to dummy atoms in the chain
-C first is connected prevous chain (itype(i+1).eq.ntyp1)=true
-C second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif !fail
- do j=1,3
- c(j,i)=c(j,i-1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i-2)-c(j,i-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i-1)+dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- else !itype(i+1).eq.ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the first dummy residue
- call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,i)=c(j,i+1)-1.9d0*e2(j)
- enddo
- else !unres_pdb
- do j=1,3
- dcj=(c(j,i+3)-c(j,i+2))/2.0
- if (dcj.eq.0) dcj=1.23591524223
- c(j,i)=c(j,i+1)-dcj
- c(j,nres+i)=c(j,i)
- enddo
- endif !unres_pdb
- endif !itype(i+1).eq.ntyp1
- endif !itype.eq.ntyp1
- enddo
+ 10 continue
+#ifdef DEBUG
+ write (iout,'(a,i5)') ' Number of residues found: ',ires
+#endif
+ if (ires.eq.0) return
C Calculate the CM of the last side chain.
+ if (iii.gt.0) then
if (unres_pdb) then
do j=1,3
dc(j,ires)=sccor(j,iii)
else
call sccenter(ires,iii,sccor)
endif
+ endif
+ nres=ires
nsup=nres
nstart_sup=1
if (itype(nres).ne.10) then
nres=nres+1
itype(nres)=ntyp1
- if (unres_pdb) then
-C 2/15/2013 by Adam: corrected insertion of the last dummy residue
- call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
- if (fail) then
- e2(1)=0.0d0
- e2(2)=1.0d0
- e2(3)=0.0d0
- endif
- do j=1,3
- c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
- enddo
- else
do j=1,3
- dcj=(c(j,nres-2)-c(j,nres-3))/2.0
- if (dcj.eq.0) dcj=1.23591524223
+ dcj=c(j,nres-2)-c(j,nres-3)
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
endif
- endif
do i=2,nres-1
do j=1,3
c(j,i+nres)=dc(j,i)
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-1.9d0*e2(j)
+ c(j,1)=c(j,2)-3.8d0*e2(j)
enddo
else
do j=1,3
- dcj=(c(j,4)-c(j,3))/2.0
+ dcj=c(j,4)-c(j,3)
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
endif
endif
-C Copy the coordinates to reference coordinates
-c do i=1,2*nres
-c do j=1,3
-c cref(j,i)=c(j,i)
-c enddo
-c enddo
C Calculate internal coordinates.
- if (out_template_coord) then
+ if (lprn) then
write (iout,'(/a)')
& "Cartesian coordinates of the reference structure"
write (iout,'(a,3(3x,a5),5x,3(3x,a5))')
enddo
endif
C Calculate internal coordinates.
-c call int_from_cart1(.false.)
- call int_from_cart(.true.,.true.)
- call sc_loc_geom(.true.)
+ write (iout,'(a)')
+ & "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)')
+ & ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),
+ & (c(j,nres+ires),j=1,3)
+ enddo
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
c write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
c & vbld_inv(i+nres)
enddo
- do i=1,nres
- do j=1,3
- cref(j,i)=c(j,i)
- cref(j,i+nres)=c(j,i+nres)
- enddo
- enddo
+c call chainbuild
+C Copy the coordinates to reference coordinates
do i=1,2*nres
do j=1,3
+ cref(j,i)=c(j,i)
chomo(j,i,k)=c(j,i)
enddo
enddo
+
+ ishift_pdb=ishift
return
end
-
-
projects / unres.git / commitdiff