From a330d50520ccfff39cfd04ca0cbd6cd85a578d6b Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Wed, 1 Apr 2020 20:03:56 +0200 Subject: [PATCH] boxshift &readpdb --- source/cluster/wham/src-HCD/DIMENSIONS | 3 + .../wham/src-HCD/Makefile-MPICH-ifort-okeanos | 4 +- source/cluster/wham/src-HCD/readpdb.F | 590 ++++++------ .../unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos | 15 +- source/unres/src-HCD-5D/boxshift.f | 101 ++ source/unres/src-HCD-5D/elecont.f | 54 +- source/unres/src-HCD-5D/energy_p_new-sep_barrier.F | 506 +++------- source/unres/src-HCD-5D/energy_p_new_barrier.F | 1016 ++++++-------------- source/unres/src-HCD-5D/gen_rand_conf.F | 19 +- source/unres/src-HCD-5D/int_from_cart.F | 269 ++++++ source/unres/src-HCD-5D/minimize_p.F | 3 +- source/unres/src-HCD-5D/readpdb-mult.F | 635 ++++++++++++ source/unres/src-HCD-5D/readpdb.F | 691 ++++--------- source/unres/src-HCD-5D/readrtns_CSA.F | 4 +- source/unres/src-HCD-5D/rmscalc.F | 4 +- source/unres/src-HCD-5D/sc_move.F | 1 + source/unres/src-HCD-5D/stochfric.F | 2 +- source/unres/src-HCD-5D/unres.F | 4 +- source/wham/src-HCD/Makefile_MPICH_ifort-okeanos | 4 +- source/wham/src-HCD/readpdb.F | 588 ++++++----- 20 files changed, 2238 insertions(+), 2275 deletions(-) create mode 100644 source/unres/src-HCD-5D/boxshift.f create mode 100644 source/unres/src-HCD-5D/int_from_cart.F create mode 100644 source/unres/src-HCD-5D/readpdb-mult.F diff --git a/source/cluster/wham/src-HCD/DIMENSIONS b/source/cluster/wham/src-HCD/DIMENSIONS index 80ac845..909354c 100644 --- a/source/cluster/wham/src-HCD/DIMENSIONS +++ b/source/cluster/wham/src-HCD/DIMENSIONS @@ -85,3 +85,6 @@ C Maximum number of SC local term fitting function coefficiants 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) diff --git a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos index 4f7f61f..d1e64a7 100644 --- a/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos +++ b/source/cluster/wham/src-HCD/Makefile-MPICH-ifort-okeanos @@ -1,7 +1,7 @@ #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 diff --git a/source/cluster/wham/src-HCD/readpdb.F b/source/cluster/wham/src-HCD/readpdb.F index dc6aa0a..2e73601 100644 --- a/source/cluster/wham/src-HCD/readpdb.F +++ b/source/cluster/wham/src-HCD/readpdb.F @@ -1,9 +1,9 @@ - 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' @@ -11,147 +11,176 @@ C geometry. 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 @@ -165,23 +194,57 @@ C Calculate the CM of the last side chain. 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) @@ -198,24 +261,30 @@ c & vbld_inv(i+nres) 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' @@ -224,56 +293,61 @@ c--------------------------------------------------------------------------- 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 @@ -283,41 +357,21 @@ c endif 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' @@ -329,7 +383,6 @@ c--------------------------------------------------------------------------- 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 @@ -338,7 +391,7 @@ c--------------------------------------------------------------------------- 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 @@ -358,7 +411,7 @@ c--------------------------------------------------------------------------- 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 @@ -372,10 +425,7 @@ c 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. @@ -399,16 +449,30 @@ c 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' @@ -420,15 +484,14 @@ c--------------------------------------------------------------------------- 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' @@ -440,7 +503,6 @@ C and convert the peptide geometry into virtual-chain geometry. 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) @@ -449,10 +511,8 @@ C and convert the peptide geometry into virtual-chain geometry. 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 @@ -463,27 +523,10 @@ c write (2,*) "UNRES_PDB",unres_pdb 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 @@ -497,7 +540,9 @@ 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)=sccor(j,iii) @@ -522,13 +567,6 @@ c write (iout,*) "BEG ires",ires 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) @@ -544,14 +582,14 @@ c write (iout,*) "New chain started",ires,ishift 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) @@ -576,60 +614,13 @@ c write (iout,*) "sidechain ",atom 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) @@ -637,31 +628,19 @@ C Calculate the CM of the last side chain. 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) @@ -683,24 +662,18 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue 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))') @@ -712,9 +685,15 @@ C Calculate internal coordinates. 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) @@ -733,19 +712,16 @@ c call int_from_cart1(.false.) 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 - - diff --git a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos index b7c9d19..9826556 100644 --- a/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos +++ b/source/unres/src-HCD-5D/Makefile_MPICH_ifort-okeanos @@ -40,7 +40,8 @@ object = unres.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \ 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 \ @@ -84,6 +85,15 @@ E0LL2Y: ${object} xdrf/libxdrf.a ${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 @@ -151,6 +161,9 @@ cartder.o : cartder.F 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 diff --git a/source/unres/src-HCD-5D/boxshift.f b/source/unres/src-HCD-5D/boxshift.f new file mode 100644 index 0000000..29d3406 --- /dev/null +++ b/source/unres/src-HCD-5D/boxshift.f @@ -0,0 +1,101 @@ + +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 diff --git a/source/unres/src-HCD-5D/elecont.f b/source/unres/src-HCD-5D/elecont.f index 690fd44..bf9056a 100644 --- a/source/unres/src-HCD-5D/elecont.f +++ b/source/unres/src-HCD-5D/elecont.f @@ -12,7 +12,7 @@ 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. @@ -53,13 +53,8 @@ c data epp / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/ 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 @@ -80,47 +75,10 @@ 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) diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F index c6c6832..96f7777 100644 --- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F @@ -85,6 +85,7 @@ c include 'COMMON.CONTACTS' 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 @@ -97,6 +98,7 @@ 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 @@ -106,9 +108,13 @@ cd & 'iend=',iend(i,iint) 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) @@ -187,6 +193,7 @@ c include 'COMMON.CONTACTS' 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 @@ -199,6 +206,7 @@ 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 @@ -210,9 +218,13 @@ cd & 'iend=',iend(i,iint) 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) @@ -285,6 +297,7 @@ C & 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 @@ -297,6 +310,7 @@ 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 @@ -304,9 +318,13 @@ c do iint=1,nint_gr(i) 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 @@ -384,6 +402,7 @@ C & 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 @@ -396,6 +415,7 @@ 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 @@ -403,9 +423,13 @@ c do iint=1,nint_gr(i) 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 @@ -478,6 +502,7 @@ C 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 @@ -499,6 +524,7 @@ 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) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -523,9 +549,13 @@ c dscj_inv=dsc_inv(itypj) 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) @@ -604,6 +634,7 @@ C 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 @@ -625,6 +656,7 @@ 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) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -649,9 +681,13 @@ c dscj_inv=dsc_inv(itypj) 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) @@ -731,6 +767,7 @@ C & 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 @@ -748,12 +785,8 @@ c do i=iatsc_s,iatsc_e 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) @@ -787,79 +820,27 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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 @@ -946,15 +927,13 @@ C 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 @@ -972,12 +951,8 @@ c do i=iatsc_s,iatsc_e 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) @@ -1011,67 +986,15 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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) @@ -1191,6 +1114,8 @@ 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) + 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) @@ -1217,9 +1142,18 @@ c dscj_inv=dsc_inv(itypj) 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) @@ -1317,6 +1251,7 @@ C & 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 @@ -1336,6 +1271,8 @@ c do i=iatsc_s,iatsc_e 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 @@ -1359,9 +1296,18 @@ c dscj_inv=dsc_inv(itypj) 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) @@ -1612,12 +1558,7 @@ C & .or. itype(i+4).eq.ntyp1 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) @@ -1641,12 +1582,7 @@ C & .or. itype(i-1).eq.ntyp1 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 @@ -1677,12 +1613,7 @@ C & .or. itype(i-1).eq.ntyp1 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) @@ -1764,6 +1695,7 @@ C------------------------------------------------------------------------------- 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 @@ -1796,44 +1728,10 @@ C print *,"WCHODZE2" 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) @@ -2734,6 +2632,7 @@ c write (iout,*) "evdwpp_short" 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" @@ -2754,12 +2653,7 @@ c do i=iatel_s_vdw,iatel_e_vdw 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) @@ -2781,43 +2675,10 @@ c do j=ielstart_vdw(i),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) @@ -2881,14 +2742,12 @@ C 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 @@ -2907,12 +2766,7 @@ c do i=iatscp_s,iatscp_e 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) @@ -2928,44 +2782,10 @@ C Uncomment following three lines for Ca-p interactions 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) @@ -3063,6 +2883,7 @@ C & 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' @@ -3076,12 +2897,7 @@ c & ' iatscp_e=',iatscp_e 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), @@ -3102,49 +2918,13 @@ C Uncomment following three lines for Ca-p interactions 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))) diff --git a/source/unres/src-HCD-5D/energy_p_new_barrier.F b/source/unres/src-HCD-5D/energy_p_new_barrier.F index ba7cbd8..ddc2a4d 100644 --- a/source/unres/src-HCD-5D/energy_p_new_barrier.F +++ b/source/unres/src-HCD-5D/energy_p_new_barrier.F @@ -115,6 +115,7 @@ c call chainbuild_cart 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 @@ -1476,6 +1477,7 @@ C & 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 @@ -1488,6 +1490,7 @@ 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 @@ -1499,9 +1502,13 @@ cd & 'iend=',iend(i,iint) 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 @@ -1649,6 +1656,7 @@ C & 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 @@ -1661,6 +1669,7 @@ 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 @@ -1668,9 +1677,13 @@ c do iint=1,nint_gr(i) 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 @@ -1748,6 +1761,7 @@ C 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 @@ -1769,6 +1783,7 @@ 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) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1803,9 +1818,13 @@ c chip12=0.0D0 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) @@ -1882,14 +1901,13 @@ C 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 @@ -1912,66 +1930,14 @@ c do i=iatsc_s,iatsc_e 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 @@ -1998,15 +1964,15 @@ c write(iout,*) "PO ZWYKLE", evdwij 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 @@ -2014,30 +1980,30 @@ 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 @@ -2048,191 +2014,112 @@ c chip12=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 @@ -2262,7 +2149,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.CALC' include 'COMMON.SPLITELE' - integer xshift,yshift,zshift,subchap + double precision boxshift integer icall common /srutu/ icall logical lprn @@ -2271,8 +2158,7 @@ C 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 @@ -2290,41 +2176,14 @@ c do i=iatsc_s,iatsc_e 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) @@ -2361,131 +2220,77 @@ c chip12=0.0D0 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 @@ -2636,6 +2441,7 @@ C 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 @@ -2648,6 +2454,7 @@ 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 @@ -2657,9 +2464,9 @@ cd & 'iend=',iend(i,iint) 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) @@ -2716,7 +2523,7 @@ c include 'COMMON.CONTACTS' 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 @@ -2732,12 +2539,7 @@ C write(iout,*) 'In EELEC_soft_sphere' 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) @@ -2754,43 +2556,10 @@ c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',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) @@ -3786,12 +3555,7 @@ c end if 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) @@ -3846,13 +3610,7 @@ c if ((zmedi.gt.((0.5d0)*boxzsize)).or. 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 @@ -3895,12 +3653,7 @@ c & .or. itype(i-1).eq.ntyp1 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 @@ -3940,7 +3693,7 @@ c do j=ielstart(i),ielend(i) 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) @@ -3948,7 +3701,7 @@ C end of changes by Ana 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 @@ -4017,10 +3770,9 @@ C------------------------------------------------------------------------------- & 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/ @@ -4032,7 +3784,6 @@ C 13-go grudnia roku pamietnego... double precision unmat(3,3) /1.0d0,0.0d0,0.0d0, & 0.0d0,1.0d0,0.0d0, & 0.0d0,0.0d0,1.0d0/ - integer xshift,yshift,zshift c time00=MPI_Wtime() cd write (iout,*) "eelecij",i,j c ind=ind+1 @@ -4055,73 +3806,12 @@ C zj=c(3,j)+0.5D0*dzj-zmedi 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) @@ -5569,7 +5259,7 @@ C 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 @@ -5612,12 +5302,7 @@ 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 xi=xi+xshift*boxxsize C yi=yi+yshift*boxysize C zi=zi+zshift*boxzsize @@ -5634,67 +5319,10 @@ C Uncomment following three lines for Ca-p interactions 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 @@ -5716,32 +5344,6 @@ C ggg(1)=xj*fac ggg(2)=yj*fac ggg(3)=zj*fac -cgrad if (j.lt.i) then -cd write (iout,*) 'ji' -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) @@ -5774,15 +5376,13 @@ C 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' @@ -5801,43 +5401,7 @@ c do i=iatscp_s,iatscp_e 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) @@ -5851,68 +5415,10 @@ C Uncomment following three lines for Ca-p interactions 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 @@ -13638,3 +13144,103 @@ C----------------------------------------------------------------------- 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 diff --git a/source/unres/src-HCD-5D/gen_rand_conf.F b/source/unres/src-HCD-5D/gen_rand_conf.F index 3e662cc..557f435 100644 --- a/source/unres/src-HCD-5D/gen_rand_conf.F +++ b/source/unres/src-HCD-5D/gen_rand_conf.F @@ -845,11 +845,12 @@ c overlapping residues left, or false otherwise (success) 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)) @@ -884,7 +885,7 @@ c write (iout,*) "i,j",i,j," itypi,itypj",itypi,itypj 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 @@ -896,17 +897,23 @@ c & ,rcomp 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 diff --git a/source/unres/src-HCD-5D/int_from_cart.F b/source/unres/src-HCD-5D/int_from_cart.F new file mode 100644 index 0000000..2ae703b --- /dev/null +++ b/source/unres/src-HCD-5D/int_from_cart.F @@ -0,0 +1,269 @@ + 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 diff --git a/source/unres/src-HCD-5D/minimize_p.F b/source/unres/src-HCD-5D/minimize_p.F index a56e4f8..cea54c4 100644 --- a/source/unres/src-HCD-5D/minimize_p.F +++ b/source/unres/src-HCD-5D/minimize_p.F @@ -635,8 +635,7 @@ c v(25)=4.0D0 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 diff --git a/source/unres/src-HCD-5D/readpdb-mult.F b/source/unres/src-HCD-5D/readpdb-mult.F new file mode 100644 index 0000000..8fd1da8 --- /dev/null +++ b/source/unres/src-HCD-5D/readpdb-mult.F @@ -0,0 +1,635 @@ + 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 + diff --git a/source/unres/src-HCD-5D/readpdb.F b/source/unres/src-HCD-5D/readpdb.F index 943d67d..9c8462d 100644 --- a/source/unres/src-HCD-5D/readpdb.F +++ b/source/unres/src-HCD-5D/readpdb.F @@ -1,7 +1,7 @@ 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' @@ -11,27 +11,33 @@ C geometry. 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. @@ -49,145 +55,112 @@ crc to be corrected !!! 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 @@ -202,12 +175,11 @@ C 2/15/2013 by Adam: corrected insertion of the last dummy residue 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 @@ -234,27 +206,46 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue 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) @@ -264,8 +255,6 @@ c write(iout,*)"before int_from_cart nres",nres 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 @@ -275,22 +264,15 @@ c & vbld_inv(i+1) 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 @@ -302,283 +284,14 @@ cc enddiag 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' @@ -588,21 +301,20 @@ C and convert the peptide geometry into virtual-chain geometry. 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 @@ -613,27 +325,10 @@ c write (2,*) "UNRES_PDB",unres_pdb 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 @@ -647,7 +342,9 @@ 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)=sccor(j,iii) @@ -672,13 +369,6 @@ c write (iout,*) "BEG ires",ires 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) @@ -694,14 +384,14 @@ c write (iout,*) "New chain started",ires,ishift 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) @@ -726,61 +416,13 @@ c write (iout,*) "sidechain ",atom 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) @@ -788,31 +430,19 @@ C Calculate the CM of the last side chain. 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) @@ -834,24 +464,18 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue 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))') @@ -863,7 +487,16 @@ C Calculate internal coordinates. 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) @@ -883,18 +516,16 @@ C Calculate internal coordinates. 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 - diff --git a/source/unres/src-HCD-5D/readrtns_CSA.F b/source/unres/src-HCD-5D/readrtns_CSA.F index 8f4d05d..0320484 100644 --- a/source/unres/src-HCD-5D/readrtns_CSA.F +++ b/source/unres/src-HCD-5D/readrtns_CSA.F @@ -177,8 +177,8 @@ c call readi(controlcard,'IZ_SC',iz_sc,0) 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) diff --git a/source/unres/src-HCD-5D/rmscalc.F b/source/unres/src-HCD-5D/rmscalc.F index b467d9c..dc976b8 100644 --- a/source/unres/src-HCD-5D/rmscalc.F +++ b/source/unres/src-HCD-5D/rmscalc.F @@ -19,6 +19,7 @@ C Loop over chain permutations 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), @@ -199,8 +200,9 @@ c------------------------------------------------------------------------ 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 diff --git a/source/unres/src-HCD-5D/sc_move.F b/source/unres/src-HCD-5D/sc_move.F index 75b7211..201e20f 100644 --- a/source/unres/src-HCD-5D/sc_move.F +++ b/source/unres/src-HCD-5D/sc_move.F @@ -156,6 +156,7 @@ c Put the original weights back to calculate the full energy mask_side=1 crc n_fun=n_fun+1 ct write (iout,*) 'sc_local time= ',MPI_WTIME()-wtime + call chainbuild_extconf return end diff --git a/source/unres/src-HCD-5D/stochfric.F b/source/unres/src-HCD-5D/stochfric.F index 09b6877..c8a3a0d 100644 --- a/source/unres/src-HCD-5D/stochfric.F +++ b/source/unres/src-HCD-5D/stochfric.F @@ -525,7 +525,7 @@ c save licznik 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 diff --git a/source/unres/src-HCD-5D/unres.F b/source/unres/src-HCD-5D/unres.F index f556eb6..08caade 100644 --- a/source/unres/src-HCD-5D/unres.F +++ b/source/unres/src-HCD-5D/unres.F @@ -283,9 +283,9 @@ crc overlap test 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 diff --git a/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos b/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos index e667382..5280f0a 100644 --- a/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos +++ b/source/wham/src-HCD/Makefile_MPICH_ifort-okeanos @@ -1,9 +1,9 @@ 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 diff --git a/source/wham/src-HCD/readpdb.F b/source/wham/src-HCD/readpdb.F index b8ce4f4..cde2738 100644 --- a/source/wham/src-HCD/readpdb.F +++ b/source/wham/src-HCD/readpdb.F @@ -1,10 +1,11 @@ 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' @@ -12,146 +13,176 @@ C geometry. 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 @@ -165,21 +196,57 @@ C Calculate the CM of the last side chain. 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) @@ -196,24 +263,30 @@ c & vbld_inv(i+nres) 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' @@ -223,56 +296,62 @@ c--------------------------------------------------------------------------- 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 @@ -282,41 +361,21 @@ c endif 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' @@ -329,7 +388,6 @@ c--------------------------------------------------------------------------- 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 @@ -338,7 +396,7 @@ c--------------------------------------------------------------------------- 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 @@ -358,7 +416,7 @@ c--------------------------------------------------------------------------- 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 @@ -372,10 +430,7 @@ c 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. @@ -399,16 +454,30 @@ c 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' @@ -420,19 +489,19 @@ c--------------------------------------------------------------------------- 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' @@ -441,7 +510,6 @@ C and convert the peptide geometry into virtual-chain geometry. 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) @@ -450,10 +518,8 @@ C and convert the peptide geometry into virtual-chain geometry. 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 @@ -464,27 +530,10 @@ c write (2,*) "UNRES_PDB",unres_pdb 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 @@ -498,7 +547,9 @@ 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)=sccor(j,iii) @@ -523,13 +574,6 @@ c write (iout,*) "BEG ires",ires 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) @@ -545,14 +589,14 @@ c write (iout,*) "New chain started",ires,ishift 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) @@ -577,60 +621,13 @@ c write (iout,*) "sidechain ",atom 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) @@ -638,31 +635,19 @@ C Calculate the CM of the last side chain. 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) @@ -684,24 +669,18 @@ C 2/15/2013 by Adam: corrected insertion of the first dummy residue 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))') @@ -713,9 +692,15 @@ C Calculate internal coordinates. 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) @@ -734,19 +719,16 @@ c call int_from_cart1(.false.) 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 - - -- 1.7.9.5