X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fenergy.f90;h=2fc27730ff6238fb5320ac2abb2e3d1ab6f61fed;hb=bbbdc8e18680625d3004f414aec255e9ca7b7353;hp=565c695e37d745e400ab28ab2d6f5a2ca9a16e86;hpb=8d2ab9ba185dbbc31bfe3d1c66d7e1c9d632463b;p=unres4.git diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 565c695..2fc2773 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -45,6 +45,7 @@ real(kind=8),dimension(:,:,:),allocatable :: gacont !(3,maxconts,maxres) integer,dimension(:),allocatable :: ishield_list integer,dimension(:,:),allocatable :: shield_list + real(kind=8),dimension(:),allocatable :: enetube,enecavtube ! ! 12/26/95 - H-bonding contacts ! common /contacts_hb/ @@ -124,6 +125,8 @@ gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, & gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,& grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres) +!-----------------------------NUCLEIC GRADIENT + real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl ! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2) real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,& gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres) @@ -225,7 +228,10 @@ real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, & Eafmforce,ethetacnstr real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6 - +! now energies for nulceic alone parameters + real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& + ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& + ecorr3_nucl #ifdef MPI real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw ! shielding effect varibles for MPI @@ -356,6 +362,7 @@ if (shield_mode.eq.2) then call set_shield_fac2 endif + print *,"AFTER EGB",ipot,evdw !mc !mc Sep-06: egb takes care of dynamic ss bonds too !mc @@ -421,6 +428,7 @@ ! Calculate the bond-stretching energy ! call ebond(estr) + print *,"EBOND",estr ! write(iout,*) "in etotal afer ebond",ipot ! @@ -537,7 +545,10 @@ else etube=0.0d0 endif - +!-------------------------------------------------------- + call ebond_nucl(estr_nucl) + call ebend_nucl(ebe_nucl) + print *,"after ebend", ebe_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -583,6 +594,21 @@ energia(23)=Eafmforce energia(24)=ethetacnstr energia(25)=etube +!--------------------------------------------------------------- + energia(26)=evdwpp + energia(27)=eespp + energia(28)=evdwpsb + energia(29)=eelpsb + energia(30)=evdwsb + energia(31)=eelsb + energia(32)=estr_nucl + energia(33)=ebe_nucl + energia(34)=esbloc + energia(35)=etors_nucl + energia(36)=etors_d_nucl + energia(37)=ecorr_nucl + energia(38)=ecorr3_nucl +!---------------------------------------------------------------------- ! Here are the energies showed per procesor if the are more processors ! per molecule then we sum it up in sum_energy subroutine ! print *," Processor",myrank," calls SUM_ENERGY" @@ -625,6 +651,10 @@ real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot, & eliptran,etube, Eafmforce,ethetacnstr + real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& + ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& + ecorr3_nucl + integer :: i #ifdef MPI integer :: ierr @@ -688,6 +718,9 @@ Eafmforce=energia(23) ethetacnstr=energia(24) etube=energia(25) + estr_nucl=energia(32) + ebe_nucl=energia(33) + #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & @@ -695,7 +728,8 @@ +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube& - +Eafmforce+ethetacnstr + +Eafmforce+ethetacnstr & + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -703,8 +737,8 @@ +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube& - +Eafmforce+ethetacnstr - + +Eafmforce+ethetacnstr & + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl #endif energia(0)=etot ! detecting NaNQ @@ -829,6 +863,9 @@ real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,& etube,ethetacnstr,Eafmforce + real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& + ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& + ecorr3_nucl etot=energia(0) evdw=energia(1) @@ -862,6 +899,9 @@ Eafmforce=energia(23) ethetacnstr=energia(24) etube=energia(25) + estr_nucl=energia(32) + ebe_nucl=energia(33) + #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -870,7 +910,9 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,& edihcnstr,ethetacnstr,ebr*nss,& - Uconst,eliptran,wliptran,Eafmforce,etube,wtube,etot + Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein + estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, & + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -898,6 +940,8 @@ 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/& 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & + 'ESTR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ & + 'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -907,7 +951,9 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,& ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, & - etube,wtube,etot + etube,wtube, & + estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,& + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -934,6 +980,8 @@ 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & + 'ESTR_nucl= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ & + 'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -974,9 +1022,9 @@ ! allocate(gacont(3,nres/4,iatsc_s:iatsc_e)) !(3,maxconts,maxres) do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) if (itypi.eq.ntyp1) cycle - itypi1=iabs(itype(i+1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -989,7 +1037,7 @@ !d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), !d & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -1006,7 +1054,7 @@ !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) !d write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') -!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj), !d & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm, !d & (c(k,i),k=1,3),(c(k,j),k=1,3) evdw=evdw+evdwij @@ -1132,9 +1180,9 @@ ! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) if (itypi.eq.ntyp1) cycle - itypi1=iabs(itype(i+1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1143,7 +1191,7 @@ ! do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -1161,7 +1209,7 @@ !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) !d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)') -!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj), !d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, !d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, !d & (c(k,i),k=1,3),(c(k,j),k=1,3) @@ -1233,9 +1281,9 @@ ! endif !el ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) if (itypi.eq.ntyp1) cycle - itypi1=iabs(itype(i+1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1250,7 +1298,7 @@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) !el ind=ind+1 - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -1302,7 +1350,7 @@ sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),15(0pf7.3))') -!d & restyp(itypi),i,restyp(itypj),j, +!d & restyp(itypi,1),i,restyp(itypj,1),j, !d & epsi,sigm,chi1,chi2,chip1,chip2, !d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq), !d & om1,om2,om12,1.0D0/dsqrt(rrij), @@ -1365,10 +1413,10 @@ !el ind=0 do i=iatsc_s,iatsc_e !C print *,"I am in EVDW",i - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) ! if (i.ne.47) cycle if (itypi.eq.ntyp1) cycle - itypi1=iabs(itype(i+1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1444,14 +1492,14 @@ enddo! k ELSE !el ind=ind+1 - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle ! if (j.ne.78) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) ! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,& ! 1.0d0/vbld(j+nres) !d -! write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -1573,7 +1621,7 @@ if (rij_shift.le.0.0D0) then evdw=1.0D20 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') -!d & restyp(itypi),i,restyp(itypj),j, +!d & restyp(itypi,1),i,restyp(itypj,1),j, !d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) return endif @@ -1596,7 +1644,7 @@ sigm=dabs(aa/bb)**(1.0D0/6.0D0) epsi=bb**2/aa!(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & - restyp(itypi),i,restyp(itypj),j, & + restyp(itypi,1),i,restyp(itypj,1),j, & epsi,sigm,chi1,chi2,chip1,chip2, & eps1,eps2rt**2,eps3rt**2,sig,sig0ij, & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift, & @@ -1608,6 +1656,7 @@ !C print *,i,j,c(1,i),c(1,j),c(2,i),c(2,j),c(3,i),c(3,j) ! if (energy_dec) write (iout,*) & ! 'evdw',i,j,evdwij +! print *,"ZALAMKA", evdw ! Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1640,6 +1689,7 @@ enddo ! j enddo ! iint enddo ! i +! print *,"ZALAMKA", evdw ! write (iout,*) "Number of loop steps in EGB:",ind !ccc energy_dec=.false. return @@ -1678,9 +1728,9 @@ ! if (icall.eq.0) lprn=.true. !el ind=0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) if (itypi.eq.ntyp1) cycle - itypi1=iabs(itype(i+1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1695,7 +1745,7 @@ do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) !el ind=ind+1 - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -1757,7 +1807,7 @@ bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & - restyp(itypi),i,restyp(itypj),j,& + restyp(itypi,1),i,restyp(itypj,1),j,& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),& chi1,chi2,chip1,chip2,& eps1,eps2rt**2,eps3rt**2,& @@ -1810,9 +1860,9 @@ evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) if (itypi.eq.ntyp1) cycle - itypi1=iabs(itype(i+1)) + itypi1=iabs(itype(i+1,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -1823,7 +1873,7 @@ !d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), !d & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -1897,7 +1947,7 @@ eello_turn4=0.0d0 !el ind=0 do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -1907,7 +1957,7 @@ num_conti=0 ! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle + if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle !el ind=ind+1 iteli=itel(i) itelj=itel(j) @@ -2340,17 +2390,17 @@ endif ! if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then if (i.gt. nnt+2 .and. i.lt.nct+2) then - iti = itortyp(itype(i-2)) + iti = itortyp(itype(i-2,1)) else iti=ntortyp+1 endif ! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - iti1 = itortyp(itype(i-1)) + iti1 = itortyp(itype(i-1,1)) else iti1=ntortyp+1 endif -! print *,iti,i,"iti",iti1,itype(i-1),itype(i-2) +! print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1) !d write (iout,*) '*******i',i,' iti1',iti !d write (iout,*) 'b1',b1(:,iti) !d write (iout,*) 'b2',b2(:,iti) @@ -2387,8 +2437,8 @@ enddo ! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then - if (itype(i-1).le.ntyp) then - iti1 = itortyp(itype(i-1)) + if (itype(i-1,1).le.ntyp) then + iti1 = itortyp(itype(i-1,1)) else iti1=ntortyp+1 endif @@ -2687,7 +2737,7 @@ #endif #endif !d do i=1,nres -!d iti = itortyp(itype(i)) +!d iti = itortyp(itype(i,1)) !d write (iout,*) i !d do j=1,2 !d write (iout,'(2f10.5,5x,2f10.5,5x,2f10.5)') @@ -2851,8 +2901,8 @@ ! print *,"before iturn3 loop" do i=iturn3_start,iturn3_end - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & - .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 & + .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2897,9 +2947,9 @@ num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & - .or. itype(i+3).eq.ntyp1 & - .or. itype(i+4).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 & + .or. itype(i+3,1).eq.ntyp1 & + .or. itype(i+4,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2940,7 +2990,7 @@ num_conti=num_cont_hb(i) call eelecij(i,i+3,ees,evdw1,eel_loc) - if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & + if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i @@ -2948,7 +2998,7 @@ ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 ! do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -2990,8 +3040,8 @@ ! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) -! write (iout,*) i,j,itype(i),itype(j) - if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle +! write (iout,*) i,j,itype(i,1),itype(j,1) + if (itype(j,1).eq.ntyp1.or. itype(j+1,1).eq.ntyp1) cycle call eelecij(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -3493,7 +3543,7 @@ a32=a32*fac a33=a33*fac !d write (iout,'(4i5,4f10.5)') -!d & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33 +!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33 !d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij !d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i), !d & uy(:,j),uz(:,j) @@ -4326,9 +4376,9 @@ a_temp(1,2)=a23 a_temp(2,1)=a32 a_temp(2,2)=a33 - iti1=itortyp(itype(i+1)) - iti2=itortyp(itype(i+2)) - iti3=itortyp(itype(i+3)) + iti1=itortyp(itype(i+1,1)) + iti2=itortyp(itype(i+2,1)) + iti3=itortyp(itype(i+3,1)) ! write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3 call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) @@ -4597,7 +4647,7 @@ !d print '(a)','Enter ESCP' !d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do i=iatscp_s,iatscp_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) @@ -4606,8 +4656,8 @@ do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - if (itype(j).eq.ntyp1) cycle - itypj=iabs(itype(j)) + if (itype(j,1).eq.ntyp1) cycle + itypj=iabs(itype(j,1)) ! Uncomment following three lines for SC-p interactions ! xj=c(1,nres+j)-xi ! yj=c(2,nres+j)-yi @@ -4701,7 +4751,7 @@ !d print '(a)','Enter ESCP' !d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do i=iatscp_s,iatscp_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) @@ -4716,7 +4766,7 @@ do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) if (itypj.eq.ntyp1) cycle ! Uncomment following three lines for SC-p interactions ! xj=c(1,nres+j)-xi @@ -4892,8 +4942,8 @@ ! 18/07/06 MC: Use the convention that the first nss pairs are SS bonds if (.not.dyn_ss .and. i.le.nss) then ! 15/02/13 CC dynamic SSbond - additional check - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. & - iabs(itype(jjj)).eq.1) then + if (ii.gt.nres .and. iabs(itype(iii,1)).eq.1 .and. & + iabs(itype(jjj,1)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij !d write (iout,*) "eij",eij @@ -5021,7 +5071,7 @@ deltat1,deltat2,deltat12,ed,pom1,pom2,eom1,eom2,eom12,& cosphi,ggk - itypi=iabs(itype(i)) + itypi=iabs(itype(i,1)) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -5030,7 +5080,7 @@ dzi=dc_norm(3,nres+i) ! dsci_inv=dsc_inv(itypi) dsci_inv=vbld_inv(nres+i) - itypj=iabs(itype(j)) + itypj=iabs(itype(j,1)) ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(nres+j) xj=c(1,nres+j)-xi @@ -5120,8 +5170,8 @@ ! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres) do i=ibondp_start,ibondp_end - if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle - if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then + if (itype(i-1,1).eq.ntyp1 .and. itype(i,1).eq.ntyp1) cycle + if (itype(i-1,1).eq.ntyp1 .or. itype(i,1).eq.ntyp1) then !C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) !C do j=1,3 !C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) & @@ -5143,11 +5193,13 @@ ! endif enddo estr=0.5d0*AKP*estr+estr1 +! print *,"estr_bb",estr,AKP ! ! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included ! do i=ibond_start,ibond_end - iti=iabs(itype(i)) + iti=iabs(itype(i,1)) + if (iti.eq.0) print *,"WARNING WRONG SETTTING",i if (iti.ne.10 .and. iti.ne.ntyp1) then nbi=nbondterm(iti) if (nbi.eq.1) then @@ -5156,6 +5208,7 @@ "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,& AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff +! print *,"estr_sc",estr do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) enddo @@ -5184,6 +5237,11 @@ usumsqder=usumsqder+ud(j)*uprod2 enddo estr=estr+uprod/usum +! print *,"estr_sc",estr,i + + if (energy_dec) write (iout,*) & + "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,& + AKSC(1,iti),uprod/usum do j=1,3 gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) enddo @@ -5233,24 +5291,24 @@ etheta=0.0D0 ! write (*,'(a,i2)') 'EBEND ICG=',icg do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle + if (itype(i-1,1).eq.ntyp1) cycle ! Zero the energy function and its derivative at 0 or pi. call splinthet(theta(i),0.5d0*delta,ss,ssd) - it=itype(i-1) - ichir1=isign(1,itype(i-2)) - ichir2=isign(1,itype(i)) - if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1)) - if (itype(i).eq.10) ichir2=isign(1,itype(i-1)) - if (itype(i-1).eq.10) then - itype1=isign(10,itype(i-2)) - ichir11=isign(1,itype(i-2)) - ichir12=isign(1,itype(i-2)) - itype2=isign(10,itype(i)) - ichir21=isign(1,itype(i)) - ichir22=isign(1,itype(i)) + it=itype(i-1,1) + ichir1=isign(1,itype(i-2,1)) + ichir2=isign(1,itype(i,1)) + if (itype(i-2,1).eq.10) ichir1=isign(1,itype(i-1,1)) + if (itype(i,1).eq.10) ichir2=isign(1,itype(i-1,1)) + if (itype(i-1,1).eq.10) then + itype1=isign(10,itype(i-2,1)) + ichir11=isign(1,itype(i-2,1)) + ichir12=isign(1,itype(i-2,1)) + itype2=isign(10,itype(i,1)) + ichir21=isign(1,itype(i,1)) + ichir22=isign(1,itype(i,1)) endif - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-2,1).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -5263,7 +5321,7 @@ y(1)=0.0D0 y(2)=0.0D0 endif - if (i.lt.nres .and. itype(i).ne.ntyp1) then + if (i.lt.nres .and. itype(i,1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -5460,27 +5518,27 @@ etheta=0.0D0 do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle - if (itype(i-2).eq.ntyp1.or.itype(i).eq.ntyp1) cycle - if (iabs(itype(i+1)).eq.20) iblock=2 - if (iabs(itype(i+1)).ne.20) iblock=1 + if (itype(i-1,1).eq.ntyp1) cycle + if (itype(i-2,1).eq.ntyp1.or.itype(i,1).eq.ntyp1) cycle + if (iabs(itype(i+1,1)).eq.20) iblock=2 + if (iabs(itype(i+1,1)).ne.20) iblock=1 dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 theti2=0.5d0*theta(i) - ityp2=ithetyp((itype(i-1))) + ityp2=ithetyp((itype(i-1,1))) do k=1,nntheterm coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then + if (i.gt.3 .and. itype(max0(i-3,1),1).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 #else phii=phi(i) #endif - ityp1=ithetyp((itype(i-2))) + ityp1=ithetyp((itype(i-2,1))) ! propagation of chirality for glycine type do k=1,nsingle cosph1(k)=dcos(k*phii) @@ -5488,13 +5546,13 @@ enddo else phii=0.0d0 - ityp1=ithetyp(itype(i-2)) + ityp1=ithetyp(itype(i-2,1)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i+1).ne.ntyp1) then + if (i.lt.nres .and. itype(i+1,1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -5502,14 +5560,14 @@ #else phii1=phi(i+1) #endif - ityp3=ithetyp((itype(i))) + ityp3=ithetyp((itype(i,1))) do k=1,nsingle cosph2(k)=dcos(k*phii1) sinph2(k)=dsin(k*phii1) enddo else phii1=0.0d0 - ityp3=ithetyp(itype(i)) + ityp3=ithetyp(itype(i,1)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -5698,7 +5756,7 @@ escloc=0.0D0 ! write (iout,'(a)') 'ESC' do i=loc_start,loc_end - it=itype(i) + it=itype(i,1) if (it.eq.ntyp1) cycle if (it.eq.10) goto 1 nlobit=nlob(iabs(it)) @@ -6028,7 +6086,7 @@ delta=0.02d0*pi escloc=0.0D0 do i=loc_start,loc_end - if (itype(i).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1) cycle 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))) @@ -6037,7 +6095,7 @@ cosfac=dsqrt(cosfac2) sinfac2=0.5d0/(1.0d0-costtab(i+1)) sinfac=dsqrt(sinfac2) - it=iabs(itype(i)) + it=iabs(itype(i,1)) if (it.eq.10) goto 1 ! ! Compute the axes of tghe local cartesian coordinates system; store in @@ -6055,7 +6113,7 @@ y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac enddo do j = 1,3 - z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i))) + z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i,1))) enddo ! write (2,*) "i",i ! write (2,*) "x_prime",(x_prime(j),j=1,3) @@ -6087,7 +6145,7 @@ ! Compute the energy of the ith side cbain ! ! write (2,*) "xx",xx," yy",yy," zz",zz - it=iabs(itype(i)) + it=iabs(itype(i,1)) do j = 1,65 x(j) = sc_parmin(j,it) enddo @@ -6095,7 +6153,7 @@ !c diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsign(1.0,dfloat(itype(i)))*dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0,dfloat(itype(i,1)))*dsin(alph(2))*dsin(omeg(2)) write(2,'(3f8.1,3f9.3,1x,3f9.3)') & alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,& xx1,yy1,zz1 @@ -6137,7 +6195,7 @@ ! & dscp1,dscp2,sumene ! sumene = enesc(x,xx,yy,zz,cost2tab(i+1),sint2tab(i+1)) escloc = escloc + sumene -! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i) +! write (2,*) "i",i," escloc",sumene,escloc,it,itype(i,1) ! & ,zz,xx,yy !#define DEBUG #ifdef DEBUG @@ -6183,7 +6241,7 @@ ! ! Compute the gradient of esc ! -! zz=zz*dsign(1.0,dfloat(itype(i))) +! zz=zz*dsign(1.0,dfloat(itype(i,1))) pom_s1=(1.0d0+x(63))/(0.1d0 + dscp1)**2 pom_s16=6*(1.0d0+x(64))/(0.1d0 + dscp1**6)**2 pom_s2=(1.0d0+x(65))/(0.1d0 + dscp2)**2 @@ -6208,7 +6266,7 @@ +(sumene2x+sumene4x*cost2tab(i+1))*(s2+s2_6) & +(pom1+pom2)*pom_dx #ifdef DEBUG - write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i) + write(2,*), "de_dxx = ", de_dxx,de_dxx_num,itype(i,1) #endif ! sumene1y=x(3) + 2*x(6)*yy + x(9)*xx + x(10)*zz @@ -6223,7 +6281,7 @@ +(sumene2y+sumene4y*cost2tab(i+1))*(s2+s2_6) & +(pom1-pom2)*pom_dy #ifdef DEBUG - write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i) + write(2,*), "de_dyy = ", de_dyy,de_dyy_num,itype(i,1) #endif ! de_dzz =(x(24) +2*x(27)*zz +x(28)*xx +x(30)*yy & @@ -6235,14 +6293,14 @@ +x(60)*xx*yy)*cost2tab(i+1)*(s2+s2_6) & + ( x(14) + 2*x(17)*zz+ x(18)*xx + x(20)*yy)*(s2+s2_6) #ifdef DEBUG - write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i) + write(2,*), "de_dzz = ", de_dzz,de_dzz_num,itype(i,1) #endif ! de_dt = 0.5d0*sumene3*cost2tab(i+1)*(s1+s1_6) & -0.5d0*sumene4*sint2tab(i+1)*(s2+s2_6) & +pom1*pom_dt1+pom2*pom_dt2 #ifdef DEBUG - write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i) + write(2,*), "de_dt = ", de_dt,de_dt_num,itype(i,1) #endif ! ! @@ -6269,9 +6327,9 @@ dZZ_Ci(k)=0.0d0 do j=1,3 dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1) & - *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) + *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres) dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1) & - *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres) + *dsign(1.0d0,dfloat(itype(i,1)))*dC_norm(j,i+nres) enddo dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres)) @@ -6471,10 +6529,10 @@ etors=0.0D0 do i=iphi_start,iphi_end etors_ii=0.0D0 - if (itype(i-2).eq.ntyp1.or. itype(i-1).eq.ntyp1 & - .or. itype(i).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 & + .or. itype(i,1).eq.ntyp1) cycle + itori=itortyp(itype(i-2,1)) + itori1=itortyp(itype(i-1,1)) phii=phi(i) gloci=0.0D0 ! Proline-Proline pair is a special case... @@ -6514,7 +6572,7 @@ 'etor',i,etors_ii if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & - restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,& + restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,& (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci ! write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) @@ -6574,17 +6632,17 @@ ! lprn=.true. etors=0.0D0 do i=iphi_start,iphi_end - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 & - .or. itype(i-3).eq.ntyp1 & - .or. itype(i).eq.ntyp1) cycle + if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 & + .or. itype(i-3,1).eq.ntyp1 & + .or. itype(i,1).eq.ntyp1) cycle etors_ii=0.0D0 - if (iabs(itype(i)).eq.20) then + if (iabs(itype(i,1)).eq.20) then iblock=2 else iblock=1 endif - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) + itori=itortyp(itype(i-2,1)) + itori1=itortyp(itype(i-1,1)) phii=phi(i) gloci=0.0D0 ! Regular cosine and sine terms @@ -6623,7 +6681,7 @@ 'etor',i,etors_ii-v0(itori,itori1,iblock) if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & - restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,& + restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,itori,itori1,& (v1(j,itori,itori1,iblock),j=1,6),& (v2(j,itori,itori1,iblock),j=1,6) gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci @@ -6685,18 +6743,18 @@ ! write(iout,*) "a tu??" do i=iphid_start,iphid_end etors_d_ii=0.0D0 - if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 & - .or. itype(i-3).eq.ntyp1 & - .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle - itori=itortyp(itype(i-2)) - itori1=itortyp(itype(i-1)) - itori2=itortyp(itype(i)) + if (itype(i-2,1).eq.ntyp1 .or. itype(i-1,1).eq.ntyp1 & + .or. itype(i-3,1).eq.ntyp1 & + .or. itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle + itori=itortyp(itype(i-2,1)) + itori1=itortyp(itype(i-1,1)) + itori2=itortyp(itype(i,1)) phii=phi(i) phii1=phi(i+1) gloci1=0.0D0 gloci2=0.0D0 iblock=1 - if (iabs(itype(i+1)).eq.20) iblock=2 + if (iabs(itype(i+1,1)).eq.20) iblock=2 ! Regular cosine and sine terms do j=1,ntermd_1(itori,itori1,itori2,iblock) @@ -6776,10 +6834,10 @@ ! write (iout,*) "EBACK_SC_COR",itau_start,itau_end esccor=0.0D0 do i=itau_start,itau_end - if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle + if ((itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1)) cycle esccor_ii=0.0D0 - isccori=isccortyp(itype(i-2)) - isccori1=isccortyp(itype(i-1)) + isccori=isccortyp(itype(i-2,1)) + isccori1=isccortyp(itype(i-1,1)) ! write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1) phii=phi(i) @@ -6791,17 +6849,17 @@ ! 2 = Ca...Ca...Ca...SC ! 3 = SC...Ca...Ca...SCi gloci=0.0D0 - if (((intertyp.eq.3).and.((itype(i-2).eq.10).or. & - (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or. & - (itype(i-1).eq.ntyp1))) & - .or. ((intertyp.eq.1).and.((itype(i-2).eq.10) & - .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1) & - .or.(itype(i).eq.ntyp1))) & - .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or. & - (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. & - (itype(i-3).eq.ntyp1)))) cycle - if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle - if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1)) & + if (((intertyp.eq.3).and.((itype(i-2,1).eq.10).or. & + (itype(i-1,1).eq.10).or.(itype(i-2,1).eq.ntyp1).or. & + (itype(i-1,1).eq.ntyp1))) & + .or. ((intertyp.eq.1).and.((itype(i-2,1).eq.10) & + .or.(itype(i-2,1).eq.ntyp1).or.(itype(i-1,1).eq.ntyp1) & + .or.(itype(i,1).eq.ntyp1))) & + .or.((intertyp.eq.2).and.((itype(i-1,1).eq.10).or. & + (itype(i-1,1).eq.ntyp1).or.(itype(i-2,1).eq.ntyp1).or. & + (itype(i-3,1).eq.ntyp1)))) cycle + if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1,1).eq.ntyp1)) cycle + if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres,1).eq.ntyp1)) & cycle do j=1,nterm_sccor(isccori,isccori1) v1ij=v1sccor(j,intertyp,isccori,isccori1) @@ -6818,7 +6876,7 @@ gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci if (lprn) & write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)') & - restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,isccori,isccori1,& + restyp(itype(i-2,1),1),i-2,restyp(itype(i-1,1),1),i-1,isccori,isccori1,& (v1sccor(j,intertyp,isccori,isccori1),j=1,6),& (v2sccor(j,intertyp,isccori,isccori1),j=1,6) gsccor_loc(i-3)=gsccor_loc(i-3)+gloci @@ -7938,9 +7996,9 @@ allocate(dipderx(3,5,4,maxconts,nres)) ! - iti1 = itortyp(itype(i+1)) + iti1 = itortyp(itype(i+1,1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1)) + itj1 = itortyp(itype(j+1,1)) else itj1=ntortyp+1 endif @@ -8033,14 +8091,14 @@ if (l.eq.j+1) then ! parallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itortyp(itype(i,1)) else iti=ntortyp+1 endif - itk1=itortyp(itype(k+1)) - itj=itortyp(itype(j)) + itk1=itortyp(itype(k+1,1)) + itj=itortyp(itype(j,1)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itortyp(itype(l+1,1)) else itl1=ntortyp+1 endif @@ -8186,15 +8244,15 @@ else ! Antiparallel orientation of the two CA-CA-CA frames. if (i.gt.1) then - iti=itortyp(itype(i)) + iti=itortyp(itype(i,1)) else iti=ntortyp+1 endif - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + itk1=itortyp(itype(k+1,1)) + itl=itortyp(itype(l,1)) + itj=itortyp(itype(j,1)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itortyp(itype(j+1,1)) else itj1=ntortyp+1 endif @@ -8542,9 +8600,9 @@ !d write (iout,*) !d & 'EELLO5: Contacts have occurred for peptide groups',i,j, !d & ' and',k,l - itk=itortyp(itype(k)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + itk=itortyp(itype(k,1)) + itl=itortyp(itype(l,1)) + itj=itortyp(itype(j,1)) eello5_1=0.0d0 eello5_2=0.0d0 eello5_3=0.0d0 @@ -9072,7 +9130,7 @@ ! i i C ! C !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - itk=itortyp(itype(k)) + itk=itortyp(itype(k,1)) s1= scalar2(AEAb1(1,2,imat),CUgb2(1,i)) s2=-scalar2(AEAb2(1,1,imat),Ug2Db1t(1,k)) s3= scalar2(AEAb2(1,1,imat),CUgb2(1,k)) @@ -9368,16 +9426,16 @@ ! ! 4/7/01 AL Component s1 was removed, because it pertains to the respective ! energy moment and not to the cluster cumulant. - iti=itortyp(itype(i)) + iti=itortyp(itype(i,1)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itortyp(itype(j+1,1)) else itj1=ntortyp+1 endif - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) + itk=itortyp(itype(k,1)) + itk1=itortyp(itype(k+1,1)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itortyp(itype(l+1,1)) else itl1=ntortyp+1 endif @@ -9489,22 +9547,22 @@ ! 4/7/01 AL Component s1 was removed, because it pertains to the respective ! energy moment and not to the cluster cumulant. !d write (2,*) 'eello_graph4: wturn6',wturn6 - iti=itortyp(itype(i)) - itj=itortyp(itype(j)) + iti=itortyp(itype(i,1)) + itj=itortyp(itype(j,1)) if (j.lt.nres-1) then - itj1=itortyp(itype(j+1)) + itj1=itortyp(itype(j+1,1)) else itj1=ntortyp+1 endif - itk=itortyp(itype(k)) + itk=itortyp(itype(k,1)) if (k.lt.nres-1) then - itk1=itortyp(itype(k+1)) + itk1=itortyp(itype(k+1,1)) else itk1=ntortyp+1 endif - itl=itortyp(itype(l)) + itl=itortyp(itype(l,1)) if (l.lt.nres-1) then - itl1=itortyp(itype(l+1)) + itl1=itortyp(itype(l+1,1)) else itl1=ntortyp+1 endif @@ -9731,11 +9789,11 @@ j=i+4 k=i+1 l=i+3 - iti=itortyp(itype(i)) - itk=itortyp(itype(k)) - itk1=itortyp(itype(k+1)) - itl=itortyp(itype(l)) - itj=itortyp(itype(j)) + iti=itortyp(itype(i,1)) + itk=itortyp(itype(k,1)) + itk1=itortyp(itype(k+1,1)) + itl=itortyp(itype(l,1)) + itj=itortyp(itype(j,1)) !d write (2,*) 'itk',itk,' itk1',itk1,' itl',itl,' itj',itj !d write (2,*) 'i',i,' k',k,' j',j,' l',l !d if (i.ne.1 .or. j.ne.3 .or. k.ne.2 .or. l.ne.4) then @@ -10280,10 +10338,8 @@ +wturn3*gshieldc_t3(j,i)& +wturn4*gshieldc_t4(j,i)& +wel_loc*gshieldc_ll(j,i)& - +wtube*gg_tube(j,i) - - - + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) enddo enddo #else @@ -10305,9 +10361,8 @@ +wcorr*gshieldc_ec(j,i) & +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i)& - +wtube*gg_tube(j,i) - - + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) enddo enddo @@ -10464,7 +10519,9 @@ +wturn4*gshieldc_loc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) & +wel_loc*gshieldc_loc_ll(j,i) & - +wtube*gg_tube(j,i) + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) + #else @@ -10498,7 +10555,9 @@ +wturn4*gshieldc_loc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) & +wel_loc*gshieldc_loc_ll(j,i) & - +wtube*gg_tube(j,i) + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) + @@ -10514,7 +10573,9 @@ +wturn3*gshieldx_t3(j,i) & +wturn4*gshieldx_t4(j,i) & +wel_loc*gshieldx_ll(j,i)& - +wtube*gg_tube_sc(j,i) + +wtube*gg_tube_sc(j,i) & + +wbond_nucl*gradbx_nucl(j,i) + enddo @@ -10724,7 +10785,7 @@ ! include 'COMMON.CALC' ! include 'COMMON.IOUNITS' real(kind=8), dimension(3) :: dcosom1,dcosom2 - +! print *,"wchodze" eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 & @@ -11100,7 +11161,7 @@ ! Derivatives in alpha and omega: ! do i=2,nres-1 -! dsci=dsc(itype(i)) +! dsci=dsc(itype(i,1)) dsci=vbld(i+nres) #ifdef OSF alphi=alph(i) @@ -12188,9 +12249,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12201,7 +12262,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), !d & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -12278,9 +12339,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12293,7 +12354,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !d write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint), !d & 'iend=',iend(i,iint) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -12368,9 +12429,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12379,7 +12440,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -12399,7 +12460,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) !d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)') -!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj), !d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, !d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, !d & (c(k,i),k=1,3),(c(k,j),k=1,3) @@ -12455,9 +12516,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12466,7 +12527,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle xj=c(1,nres+j)-xi yj=c(2,nres+j)-yi @@ -12486,7 +12547,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !d sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) !d epsi=bb(itypi,itypj)**2/aa(itypi,itypj) !d write (iout,'(2(a3,i3,2x),8(1pd12.4)/2(3(1pd12.4),5x)/)') -!d & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj), +!d & restyp(itypi,1),i,restyp(itypj,1),j,aa(itypi,itypj), !d & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm, !d & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij, !d & (c(k,i),k=1,3),(c(k,j),k=1,3) @@ -12554,9 +12615,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! endif !el ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12571,7 +12632,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) !el ind=ind+1 - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -12612,7 +12673,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),15(0pf7.3))') -!d & restyp(itypi),i,restyp(itypj),j, +!d & restyp(itypi,1),i,restyp(itypj,1),j, !d & epsi,sigm,chi1,chi2,chip1,chip2, !d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq), !d & om1,om2,om12,1.0D0/dsqrt(rrij), @@ -12674,9 +12735,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! endif !el ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12691,7 +12752,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) !el ind=ind+1 - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -12732,7 +12793,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) !d write (iout,'(2(a3,i3,2x),15(0pf7.3))') -!d & restyp(itypi),i,restyp(itypj),j, +!d & restyp(itypi,1),i,restyp(itypj,1),j, !d & epsi,sigm,chi1,chi2,chip1,chip2, !d & eps1,eps2rt**2,eps3rt**2,1.0D0/dsqrt(sigsq), !d & om1,om2,om12,1.0D0/dsqrt(rrij), @@ -12794,9 +12855,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! if (icall.eq.0) lprn=.false. !el ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -12872,13 +12933,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ELSE !el ind=ind+1 - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) ! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, ! & 1.0d0/vbld(j+nres) -! write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -12982,7 +13043,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (rij_shift.le.0.0D0) then evdw=1.0D20 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') -!d & restyp(itypi),i,restyp(itypj),j, +!d & restyp(itypi,1),i,restyp(itypj,1),j, !d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) return endif @@ -13003,7 +13064,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & - restyp(itypi),i,restyp(itypj),j,& + restyp(itypi,1),i,restyp(itypj,1),j,& epsi,sigm,chi1,chi2,chip1,chip2,& eps1,eps2rt**2,eps3rt**2,sig,sig0ij,& om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,& @@ -13074,9 +13135,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! if (icall.eq.0) lprn=.false. !el ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -13158,13 +13219,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! 'evdw',i,j,evdwij,' ss' ELSE !el ind=ind+1 - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) ! write (iout,*) "j",j,dsc_inv(itypj),dscj_inv, ! & 1.0d0/vbld(j+nres) -! write (iout,*) "i",i," j", j," itype",itype(i),itype(j) +! write (iout,*) "i",i," j", j," itype",itype(i,1),itype(j,1) sig0ij=sigma(itypi,itypj) chi1=chi(itypi,itypj) chi2=chi(itypj,itypi) @@ -13273,7 +13334,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (rij_shift.le.0.0D0) then evdw=1.0D20 !d write (iout,'(2(a3,i3,2x),17(0pf7.3))') -!d & restyp(itypi),i,restyp(itypj),j, +!d & restyp(itypi,1),i,restyp(itypj,1),j, !d & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) return endif @@ -13294,7 +13355,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & - restyp(itypi),i,restyp(itypj),j,& + restyp(itypi,1),i,restyp(itypj,1),j,& epsi,sigm,chi1,chi2,chip1,chip2,& eps1,eps2rt**2,eps3rt**2,sig,sig0ij,& om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,& @@ -13364,9 +13425,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! if (icall.eq.0) lprn=.true. !el ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -13381,7 +13442,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) !el ind=ind+1 - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -13437,7 +13498,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & - restyp(itypi),i,restyp(itypj),j,& + restyp(itypi,1),i,restyp(itypj,1),j,& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),& chi1,chi2,chip1,chip2,& eps1,eps2rt**2,eps3rt**2,& @@ -13493,9 +13554,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! if (icall.eq.0) lprn=.true. !el ind=0 do i=iatsc_s,iatsc_e - itypi=itype(i) + itypi=itype(i,1) if (itypi.eq.ntyp1) cycle - itypi1=itype(i+1) + itypi1=itype(i+1,1) xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) @@ -13510,7 +13571,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) !el ind=ind+1 - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! dscj_inv=dsc_inv(itypj) dscj_inv=vbld_inv(j+nres) @@ -13566,7 +13627,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sigm=dabs(aa_aq(itypi,itypj)/bb_aq(itypi,itypj))**(1.0D0/6.0D0) epsi=bb_aq(itypi,itypj)**2/aa_aq(itypi,itypj) write (iout,'(2(a3,i3,2x),17(0pf7.3))') & - restyp(itypi),i,restyp(itypj),j,& + restyp(itypi,1),i,restyp(itypj,1),j,& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),& chi1,chi2,chip1,chip2,& eps1,eps2rt**2,eps3rt**2,& @@ -13717,8 +13778,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Loop over i,i+2 and i,i+3 pairs of the peptide groups ! do i=iturn3_start,iturn3_end - if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1 & - .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1 & + .or. itype(i+2,1).eq.ntyp1 .or. itype(i+3,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -13740,9 +13801,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 & - .or. itype(i+3).eq.ntyp1 & - .or. itype(i+4).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1 & + .or. itype(i+3,1).eq.ntyp1 & + .or. itype(i+4,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -13760,7 +13821,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' if (zmedi.lt.0) zmedi=zmedi+boxzsize num_conti=num_cont_hb(i) call eelecij_scale(i,i+3,ees,evdw1,eel_loc) - if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) & + if (wturn4.gt.0.0d0 .and. itype(i+2,1).ne.ntyp1) & call eturn4(i,eello_turn4) num_cont_hb(i)=num_conti enddo ! i @@ -13768,7 +13829,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 ! do i=iatel_s,iatel_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -13787,7 +13848,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle + if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle call eelecij_scale(i,j,ees,evdw1,eel_loc) enddo ! j num_cont_hb(i)=num_conti @@ -14171,7 +14232,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' a32=a32*fac a33=a33*fac !d write (iout,'(4i5,4f10.5)') -!d & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33 +!d & i,itortyp(itype(i,1)),j,itortyp(itype(j,1)),a22,a23,a32,a33 !d write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij !d write (iout,'(2(3f10.5,5x)/2(3f10.5,5x))') uy(:,i),uz(:,i), !d & uy(:,j),uz(:,j) @@ -14648,7 +14709,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! & " iatel_e_vdw",iatel_e_vdw call flush(iout) do i=iatel_s_vdw,iatel_e_vdw - if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1.or. itype(i+1,1).eq.ntyp1) cycle dxi=dc(1,i) dyi=dc(2,i) dzi=dc(3,i) @@ -14669,7 +14730,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! & ' ielend',ielend_vdw(i) call flush(iout) do j=ielstart_vdw(i),ielend_vdw(i) - if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle + if (itype(j,1).eq.ntyp1 .or. itype(j+1,1).eq.ntyp1) cycle !el ind=ind+1 iteli=itel(i) itelj=itel(j) @@ -14802,7 +14863,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !d print '(a)','Enter ESCP' !d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do i=iatscp_s,iatscp_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) @@ -14817,7 +14878,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! Uncomment following three lines for SC-p interactions ! xj=c(1,nres+j)-xi @@ -14961,7 +15022,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !d print '(a)','Enter ESCP' !d write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e do i=iatscp_s,iatscp_e - if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) cycle iteli=itel(i) xi=0.5D0*(c(1,i)+c(1,i+1)) yi=0.5D0*(c(2,i)+c(2,i+1)) @@ -14976,7 +15037,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do iint=1,nscp_gr(i) do j=iscpstart(i,iint),iscpend(i,iint) - itypj=itype(j) + itypj=itype(j,1) if (itypj.eq.ntyp1) cycle ! Uncomment following three lines for SC-p interactions ! xj=c(1,nres+j)-xi @@ -15800,7 +15861,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo if (n.le.nphi+ntheta) goto 10 do i=2,nres-1 - if (itype(i).ne.10) then + if (itype(i,1).ne.10) then galphai=0.0D0 gomegai=0.0D0 do k=1,3 @@ -16096,6 +16157,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gg_tube(j,i)=0.0d0 gg_tube_sc(j,i)=0.0d0 gradafm(j,i)=0.0d0 + gradb_nucl(j,i)=0.0d0 + gradbx_nucl(j,i)=0.0d0 do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo @@ -16229,10 +16292,10 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do j=1,3 dcostheta(j,1,i)=-(dc_norm(j,i-1)+cost*dc_norm(j,i-2))/& vbld(i-1) - if (itype(i-1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint + if (itype(i-1,1).ne.ntyp1) dtheta(j,1,i)=-dcostheta(j,1,i)/sint dcostheta(j,2,i)=-(dc_norm(j,i-2)+cost*dc_norm(j,i-1))/& vbld(i) - if (itype(i-1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint + if (itype(i-1,1).ne.ntyp1) dtheta(j,2,i)=-dcostheta(j,2,i)/sint enddo enddo #if defined(MPI) && defined(PARINTDER) @@ -16241,7 +16304,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #else do i=3,nres #endif - if ((itype(i-1).ne.10).and.(itype(i-1).ne.ntyp1)) then + if ((itype(i-1,1).ne.10).and.(itype(i-1,1).ne.ntyp1)) then cost1=dcos(omicron(1,i)) sint1=sqrt(1-cost1*cost1) cost2=dcos(omicron(2,i)) @@ -16265,7 +16328,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dcosomicron(j,2,2,i)=-(dc_norm(j,i-1) & +cost2*(-dc_norm(j,i-1+nres)))/ & vbld(i-1+nres) -! write(iout,*) "vbld", i,itype(i),vbld(i-1+nres) +! write(iout,*) "vbld", i,itype(i,1),vbld(i-1+nres) domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i) enddo endif @@ -16280,7 +16343,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #else do i=4,nres #endif -! if (itype(i-1).eq.21 .or. itype(i-2).eq.21 ) cycle +! if (itype(i-1,1).eq.21 .or. itype(i-2,1).eq.21 ) cycle ! the conventional case sint=dsin(theta(i)) sint1=dsin(theta(i-1)) @@ -16305,7 +16368,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ctgt=cost/sint ctgt1=cost1/sint1 cosg_inv=1.0d0/cosg - if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then + if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then dsinphi(j,1,i)=-sing*ctgt1*dtheta(j,1,i-1) & -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2) dphi(j,1,i)=cosg_inv*dsinphi(j,1,i) @@ -16323,7 +16386,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! Obtaining the gamma derivatives from cosine derivative else do j=1,3 - if (itype(i-1).ne.ntyp1 .and. itype(i-2).ne.ntyp1) then + if (itype(i-1,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1) then dcosphi(j,1,i)=fac1*dcostheta(j,1,i-1)+fac3* & dcostheta(j,1,i-1)-fac0*(dc_norm(j,i-1)-scalp* & dc_norm(j,i-3))/vbld(i-2) @@ -16347,9 +16410,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=3,nres !elwrite(iout,*) " vecpr",i,nres #endif - if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle -! if ((itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10).or. -! & (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle + if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle +! if ((itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10).or. +! & (itype(i-1,1).eq.ntyp1).or.(itype(i,1).eq.ntyp1)) cycle !c dtauangle(j,intertyp,dervityp,residue number) !c INTERTYP=1 SC...Ca...Ca..Ca ! the conventional case @@ -16427,8 +16490,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #else do i=4,nres #endif - if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. & - (itype(i-2).eq.ntyp1).or.(itype(i-3).eq.ntyp1)) cycle + if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. & + (itype(i-2,1).eq.ntyp1).or.(itype(i-3,1).eq.ntyp1)) cycle ! the conventional case sint=dsin(omicron(1,i)) sint1=dsin(theta(i-1)) @@ -16502,8 +16565,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=3,nres #endif ! the conventional case - if ((itype(i-1).eq.ntyp1).or.(itype(i-1).eq.10).or. & - (itype(i-2).eq.ntyp1).or.(itype(i-2).eq.10)) cycle + if ((itype(i-1,1).eq.ntyp1).or.(itype(i-1,1).eq.10).or. & + (itype(i-2,1).eq.ntyp1).or.(itype(i-2,1).eq.10)) cycle sint=dsin(omicron(1,i)) sint1=dsin(omicron(2,i-1)) sing=dsin(tauangle(3,i)) @@ -16573,7 +16636,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' #else do i=2,nres-1 #endif - if(itype(i).ne.10 .and. itype(i).ne.ntyp1) then + if(itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then fac5=1.0d0/dsqrt(2*(1+dcos(theta(i+1)))) fac6=fac5/vbld(i) fac7=fac5*fac5 @@ -16807,7 +16870,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' write (iout,*) & "Analytical (upper) and numerical (lower) gradient of alpha" do i=2,nres-1 - if(itype(i).ne.10) then + if(itype(i,1).ne.10) then do j=1,3 dcji=dc(j,i-1) dc(j,i-1)=dcji+aincr @@ -16843,7 +16906,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' write (iout,*) & "Analytical (upper) and numerical (lower) gradient of omega" do i=2,nres-1 - if(itype(i).ne.10) then + if(itype(i,1).ne.10) then do j=1,3 dcji=dc(j,i-1) dc(j,i-1)=dcji+aincr @@ -16909,7 +16972,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' (cref(3,jl,kkk)-cref(3,il,kkk))**2) dij=dist(il,jl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) - if (itype(il).ne.10 .or. itype(jl).ne.10) then + if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -16936,7 +16999,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' (cref(3,jl,kkk)-cref(3,il,kkk))**2) dij=dist(il,jl) qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) - if (itype(il).ne.10 .or. itype(jl).ne.10) then + if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -16998,7 +17061,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dqwol(k,jl)=dqwol(k,jl)-ddqij enddo - if (itype(il).ne.10 .or. itype(jl).ne.10) then + if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -17039,7 +17102,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dqwol(k,il)=dqwol(k,il)+ddqij dqwol(k,jl)=dqwol(k,jl)-ddqij enddo - if (itype(il).ne.10 .or. itype(jl).ne.10) then + if (itype(il,1).ne.10 .or. itype(jl,1).ne.10) then nl=nl+1 d0ijCM=dsqrt( & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+ & @@ -17530,13 +17593,13 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !el allocate(dyn_ssbond_ij(iatsc_s:iatsc_e,nres)) !el allocate(dyn_ssbond_ij(0:nres+4,nres)) - itypi=itype(i) + itypi=itype(i,1) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) dsci_inv=vbld_inv(i+nres) - itypj=itype(j) + itypj=itype(j,1) xj=c(1,nres+j)-c(1,nres+i) yj=c(2,nres+j)-c(2,nres+i) zj=c(3,nres+j)-c(3,nres+i) @@ -17898,7 +17961,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' j=resj k=resk !C write(iout,*) resi,resj,resk - itypi=itype(i) + itypi=itype(i,1) dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -17906,7 +17969,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' xi=c(1,nres+i) yi=c(2,nres+i) zi=c(3,nres+i) - itypj=itype(j) + itypj=itype(j,1) xj=c(1,nres+j) yj=c(2,nres+j) zj=c(3,nres+j) @@ -17915,7 +17978,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' dyj=dc_norm(2,nres+j) dzj=dc_norm(3,nres+j) dscj_inv=vbld_inv(j+nres) - itypk=itype(k) + itypk=itype(k,1) xk=c(1,nres+k) yk=c(2,nres+k) zk=c(3,nres+k) @@ -18224,7 +18287,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ! print *, "I am in eliptran" do i=ilip_start,ilip_end !C do i=1,1 - if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1).or.(i.eq.nres))& + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1).or.(i.eq.nres))& cycle positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize)) @@ -18270,7 +18333,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo ! here starts the side chain transfer do i=ilip_start,ilip_end - if (itype(i).eq.ntyp1) cycle + if (itype(i,1).eq.ntyp1) cycle positi=(mod(c(3,i+nres),boxzsize)) if (positi.le.0) positi=positi+boxzsize !C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop @@ -18286,25 +18349,25 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C lipbufthick is thickenes of lipid buffore sslip=sscalelip(fracinbuf) ssgradlip=-sscagradlip(fracinbuf)/lipbufthick - eliptran=eliptran+sslip*liptranene(itype(i)) + eliptran=eliptran+sslip*liptranene(itype(i,1)) gliptranx(3,i)=gliptranx(3,i) & - +ssgradlip*liptranene(itype(i)) + +ssgradlip*liptranene(itype(i,1)) gliptranc(3,i-1)= gliptranc(3,i-1) & - +ssgradlip*liptranene(itype(i)) + +ssgradlip*liptranene(itype(i,1)) !C print *,"doing sccale for lower part" elseif (positi.gt.bufliptop) then fracinbuf=1.0d0- & ((bordliptop-positi)/lipbufthick) sslip=sscalelip(fracinbuf) ssgradlip=sscagradlip(fracinbuf)/lipbufthick - eliptran=eliptran+sslip*liptranene(itype(i)) + eliptran=eliptran+sslip*liptranene(itype(i,1)) gliptranx(3,i)=gliptranx(3,i) & - +ssgradlip*liptranene(itype(i)) + +ssgradlip*liptranene(itype(i,1)) gliptranc(3,i-1)= gliptranc(3,i-1) & - +ssgradlip*liptranene(itype(i)) + +ssgradlip*liptranene(itype(i,1)) !C print *, "doing sscalefor top part",sslip,fracinbuf else - eliptran=eliptran+liptranene(itype(i)) + eliptran=eliptran+liptranene(itype(i,1)) !C print *,"I am in true lipid" endif endif ! if in lipid or buffor @@ -18327,7 +18390,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a !C simple Kihara potential subroutine calctube(Etube) - real(kind=8) :: vectube(3),enetube(nres*2) + real(kind=8),dimension(3) :: vectube real(kind=8) :: Etube,xtemp,xminact,yminact,& ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, & sc_aa_tube,sc_bb_tube @@ -18341,7 +18404,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C for UNRES do i=itube_start,itube_end !C lets ommit dummy atoms for now - if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle !C now calculate distance from center of tube and direction vectors xmin=boxxsize ymin=boxysize @@ -18404,7 +18467,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=itube_start,itube_end !C Lets not jump over memory as we use many times iti - iti=itype(i) + iti=itype(i,1) !C lets ommit dummy atoms for now if ((iti.eq.ntyp1) & !C in UNRES uncomment the line below as GLY has no side-chain... @@ -18488,7 +18551,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C and r0 is the excluded size of nanotube (can be set to 0 if we want just a !C simple Kihara potential subroutine calctube2(Etube) - real(kind=8) :: vectube(3),enetube(nres*2) + real(kind=8),dimension(3) :: vectube real(kind=8) :: Etube,xtemp,xminact,yminact,& ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,& sstube,ssgradtube,sc_aa_tube,sc_bb_tube @@ -18504,7 +18567,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=itube_start,itube_end !C lets ommit dummy atoms for now - if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle !C now calculate distance from center of tube and direction vectors !C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) !C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize @@ -18565,12 +18628,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C lipbufthick is thickenes of lipid buffore sstube=sscalelip(fracinbuf) ssgradtube=-sscagradlip(fracinbuf)/tubebufthick -!C print *,ssgradtube, sstube,tubetranene(itype(i)) +!C print *,ssgradtube, sstube,tubetranene(itype(i,1)) enetube(i)=enetube(i)+sstube*tubetranenepep !C gg_tube_SC(3,i)=gg_tube_SC(3,i) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C gg_tube(3,i-1)= gg_tube(3,i-1) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C print *,"doing sccale for lower part" elseif (positi.gt.buftubetop) then fracinbuf=1.0d0- & @@ -18579,9 +18642,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' ssgradtube=sscagradlip(fracinbuf)/tubebufthick enetube(i)=enetube(i)+sstube*tubetranenepep !C gg_tube_SC(3,i)=gg_tube_SC(3,i) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C gg_tube(3,i-1)= gg_tube(3,i-1) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C print *, "doing sscalefor top part",sslip,fracinbuf else sstube=1.0d0 @@ -18622,7 +18685,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C print *,gg_tube(1,0),"TU" do i=itube_start,itube_end !C Lets not jump over memory as we use many times iti - iti=itype(i) + iti=itype(i,1) !C lets ommit dummy atoms for now if ((iti.eq.ntyp1) & !!C in UNRES uncomment the line below as GLY has no side-chain... @@ -18654,12 +18717,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C lipbufthick is thickenes of lipid buffore sstube=sscalelip(fracinbuf) ssgradtube=-sscagradlip(fracinbuf)/tubebufthick -!C print *,ssgradtube, sstube,tubetranene(itype(i)) - enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +!C print *,ssgradtube, sstube,tubetranene(itype(i,1)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1)) !C gg_tube_SC(3,i)=gg_tube_SC(3,i) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C gg_tube(3,i-1)= gg_tube(3,i-1) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C print *,"doing sccale for lower part" elseif (positi.gt.buftubetop) then fracinbuf=1.0d0- & @@ -18667,16 +18730,16 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' sstube=sscalelip(fracinbuf) ssgradtube=sscagradlip(fracinbuf)/tubebufthick - enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1)) !C gg_tube_SC(3,i)=gg_tube_SC(3,i) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C gg_tube(3,i-1)= gg_tube(3,i-1) -!C &+ssgradtube*tubetranene(itype(i)) +!C &+ssgradtube*tubetranene(itype(i,1)) !C print *, "doing sscalefor top part",sslip,fracinbuf else sstube=1.0d0 ssgradtube=0.0d0 - enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i,1)) !C print *,"I am in true lipid" endif else @@ -18725,12 +18788,12 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' end subroutine calctube2 !===================================================================================================================================== subroutine calcnano(Etube) - real(kind=8) :: vectube(3),enetube(nres*2), & - enecavtube(nres*2) + real(kind=8),dimension(3) :: vectube + real(kind=8) :: Etube,xtemp,xminact,yminact,& ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,& sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact - integer:: i,j,iti + integer:: i,j,iti,r Etube=0.0d0 ! print *,itube_start,itube_end,"poczatek" @@ -18743,7 +18806,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C for UNRES do i=itube_start,itube_end !C lets ommit dummy atoms for now - if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle + if ((itype(i,1).eq.ntyp1).or.(itype(i+1,1).eq.ntyp1)) cycle !C now calculate distance from center of tube and direction vectors xmin=boxxsize ymin=boxysize @@ -18826,7 +18889,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' !C fac=fac+faccav !C 667 continue endif - + if (energy_dec) write(iout,*),i,rdiff,enetube(i),enecavtube(i) do j=1,3 gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 @@ -18836,7 +18899,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' do i=itube_start,itube_end enecavtube(i)=0.0d0 !C Lets not jump over memory as we use many times iti - iti=itype(i) + iti=itype(i,1) !C lets ommit dummy atoms for now if ((iti.eq.ntyp1) & !C in UNRES uncomment the line below as GLY has no side-chain... @@ -18929,6 +18992,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac enddo + if (energy_dec) write(iout,*),i,rdiff,enetube(i+nres),enecavtube(i+nres) enddo @@ -18937,6 +19001,23 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) & +enecavtube(i+nres) enddo +! do i=1,20 +! print *,"begin", i,"a" +! do r=1,10000 +! rdiff=r/100.0d0 +! rdiff6=rdiff**6.0d0 +! sc_aa_tube=sc_aa_tube_par(i) +! sc_bb_tube=sc_bb_tube_par(i) +! enetube(i)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 +! denominator=(1.0d0+dcavtub(i)*rdiff6*rdiff6) +! enecavtube(i)= & +! (bcavtub(i)*rdiff+acavtub(i)*dsqrt(rdiff)+ccavtub(i)) & +! /denominator + +! print '(5(f10.3,1x))',rdiff,enetube(i),enecavtube(i),enecavtube(i)+enetube(i) +! enddo +! print *,"end",i,"a" +! enddo !C print *,"ETUBE", etube return end subroutine calcnano @@ -18967,14 +19048,14 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo do i=ivec_start,ivec_end !C do i=1,nres-1 -!C if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle +!C if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle ishield_list(i)=0 - if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle + if ((itype(i,1).eq.ntyp1).and.itype(i+1,1).eq.ntyp1) cycle !Cif there two consequtive dummy atoms there is no peptide group between them !C the line below has to be changed for FGPROC>1 VolumeTotal=0.0 do k=1,nres - if ((itype(k).eq.ntyp1).or.(itype(k).eq.10)) cycle + if ((itype(k,1).eq.ntyp1).or.(itype(k,1).eq.10)) cycle dist_pep_side=0.0 dist_side_calf=0.0 do j=1,3 @@ -19029,8 +19110,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo endif !C this is what is now we have the distance scaling now volume... - short=short_r_sidechain(itype(k)) - long=long_r_sidechain(itype(k)) + short=short_r_sidechain(itype(k,1)) + long=long_r_sidechain(itype(k,1)) costhet=1.0d0/dsqrt(1.0d0+short**2/dist_pep_side**2) sinthet=short/dist_pep_side*costhet !C now costhet_grad @@ -19116,7 +19197,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo fac_shield(i)=VolumeTotal*wshield+(1.0d0-wshield) -!C write(2,*) "TOTAL VOLUME",i,itype(i),fac_shield(i) +!C write(2,*) "TOTAL VOLUME",i,itype(i,1),fac_shield(i) enddo return end subroutine set_shield_fac2 @@ -19424,6 +19505,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(gg_tube_sc(3,-1:nres)) allocate(gg_tube(3,-1:nres)) allocate(gradafm(3,-1:nres)) + allocate(gradb_nucl(3,-1:nres)) + allocate(gradbx_nucl(3,-1:nres)) !(3,maxres) allocate(grad_shield_side(3,50,nres)) allocate(grad_shield_loc(3,50,nres)) @@ -19552,6 +19635,9 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(shield_list(50,nres)) allocate(dyn_ss_mask(nres)) allocate(fac_shield(nres)) + allocate(enetube(nres*2)) + allocate(enecavtube(nres*2)) + !(maxres) dyn_ss_mask(:)=.false. !---------------------- @@ -19598,6 +19684,283 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' return end subroutine alloc_ener_arrays +!----------------------------------------------------------------- + subroutine ebond_nucl(estr_nucl) +!c +!c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds +!c + + real(kind=8),dimension(3) :: u,ud + real(kind=8) :: usum,uprod,uprod1,uprod2,usumsqder + real(kind=8) :: estr_nucl,diff + integer :: iti,i,j,k,nbi + estr_nucl=0.0d0 +!C print *,"I enter ebond" + if (energy_dec) & + write (iout,*) "ibondp_start,ibondp_end",& + ibondp_nucl_start,ibondp_nucl_end + do i=ibondp_nucl_start,ibondp_nucl_end + if (itype(i-1,2).eq.ntyp1_molec(2) .or. & + itype(i,2).eq.ntyp1_molec(2)) cycle +! estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) +! do j=1,3 +! gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) +! & *dc(j,i-1)/vbld(i) +! enddo +! if (energy_dec) write(iout,*) +! & "estr1",i,vbld(i),distchainmax, +! & gnmr1(vbld(i),-1.0d0,distchainmax) + + diff = vbld(i)-vbldp0_nucl + if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),& + vbldp0_nucl,diff,AKP_nucl*diff*diff + estr_nucl=estr_nucl+diff*diff + print *,estr_nucl + do j=1,3 + gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i) + enddo +!c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) + enddo + estr_nucl=0.5d0*AKP_nucl*estr_nucl + print *,"partial sum", estr_nucl,AKP_nucl + + if (energy_dec) & + write (iout,*) "ibondp_start,ibondp_end",& + ibond_nucl_start,ibond_nucl_end + + do i=ibond_nucl_start,ibond_nucl_end +!C print *, "I am stuck",i + iti=itype(i,2) + if (iti.eq.ntyp1_molec(2)) cycle + nbi=nbondterm_nucl(iti) +!C print *,iti,nbi + if (nbi.eq.1) then + diff=vbld(i+nres)-vbldsc0_nucl(1,iti) + + if (energy_dec) & + write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, & + AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff + estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff + print *,estr_nucl + do j=1,3 + gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) + enddo + else + do j=1,nbi + diff=vbld(i+nres)-vbldsc0_nucl(j,iti) + ud(j)=aksc_nucl(j,iti)*diff + u(j)=abond0_nucl(j,iti)+0.5d0*ud(j)*diff + enddo + uprod=u(1) + do j=2,nbi + uprod=uprod*u(j) + enddo + usum=0.0d0 + usumsqder=0.0d0 + do j=1,nbi + uprod1=1.0d0 + uprod2=1.0d0 + do k=1,nbi + if (k.ne.j) then + uprod1=uprod1*u(k) + uprod2=uprod2*u(k)*u(k) + endif + enddo + usum=usum+uprod1 + usumsqder=usumsqder+ud(j)*uprod2 + enddo + estr_nucl=estr_nucl+uprod/usum + do j=1,3 + gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) + enddo + endif + enddo +!C print *,"I am about to leave ebond" + return + end subroutine ebond_nucl + +!----------------------------------------------------------------------------- + subroutine ebend_nucl(etheta_nucl) + real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm + real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle + real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble + logical :: lprn=.true., lprn1=.false. +!el local variables + integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m + real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai + real(kind=8) :: aux,etheta_nucl,ccl,ssl,scl,csl,ethetacnstr +! local variables for constrains + real(kind=8) :: difi,thetiii + integer itheta + etheta_nucl=0.0D0 + print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres + do i=ithet_nucl_start,ithet_nucl_end + if ((itype(i-1,2).eq.ntyp1_molec(2)).or.& + (itype(i-2,2).eq.ntyp1_molec(2)).or. & + (itype(i,2).eq.ntyp1_molec(2))) cycle + dethetai=0.0d0 + dephii=0.0d0 + dephii1=0.0d0 + theti2=0.5d0*theta(i) + ityp2=ithetyp_nucl(itype(i-1,2)) + do k=1,nntheterm_nucl + coskt(k)=dcos(k*theti2) + sinkt(k)=dsin(k*theti2) + enddo + if (i.gt.3 .and. itype(i-2,2).ne.ntyp1_molec(2)) then +#ifdef OSF + phii=phi(i) + if (phii.ne.phii) phii=150.0 +#else + phii=phi(i) +#endif + ityp1=ithetyp_nucl(itype(i-2,2)) + do k=1,nsingle_nucl + cosph1(k)=dcos(k*phii) + sinph1(k)=dsin(k*phii) + enddo + else + phii=0.0d0 + ityp1=nthetyp_nucl+1 + do k=1,nsingle_nucl + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + endif + + if (i.lt.nres .and. itype(i,2).ne.ntyp1_molec(2)) then +#ifdef OSF + phii1=phi(i+1) + if (phii1.ne.phii1) phii1=150.0 + phii1=pinorm(phii1) +#else + phii1=phi(i+1) +#endif + ityp3=ithetyp_nucl(itype(i,2)) + do k=1,nsingle_nucl + cosph2(k)=dcos(k*phii1) + sinph2(k)=dsin(k*phii1) + enddo + else + phii1=0.0d0 + ityp3=nthetyp_nucl+1 + do k=1,nsingle_nucl + cosph2(k)=0.0d0 + sinph2(k)=0.0d0 + enddo + endif + ethetai=aa0thet_nucl(ityp1,ityp2,ityp3) + do k=1,ndouble_nucl + do l=1,k-1 + ccl=cosph1(l)*cosph2(k-l) + ssl=sinph1(l)*sinph2(k-l) + scl=sinph1(l)*cosph2(k-l) + csl=cosph1(l)*sinph2(k-l) + cosph1ph2(l,k)=ccl-ssl + cosph1ph2(k,l)=ccl+ssl + sinph1ph2(l,k)=scl+csl + sinph1ph2(k,l)=scl-csl + enddo + enddo + if (lprn) then + write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,& + " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1 + write (iout,*) "coskt and sinkt",nntheterm_nucl + do k=1,nntheterm_nucl + write (iout,*) k,coskt(k),sinkt(k) + enddo + endif + do k=1,ntheterm_nucl + ethetai=ethetai+aathet_nucl(k,ityp1,ityp2,ityp3)*sinkt(k) + dethetai=dethetai+0.5d0*k*aathet_nucl(k,ityp1,ityp2,ityp3)& + *coskt(k) + if (lprn)& + write (iout,*) "k",k," aathet",aathet_nucl(k,ityp1,ityp2,ityp3),& + " ethetai",ethetai + enddo + if (lprn) then + write (iout,*) "cosph and sinph" + do k=1,nsingle_nucl + write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k) + enddo + write (iout,*) "cosph1ph2 and sinph2ph2" + do k=2,ndouble_nucl + do l=1,k-1 + write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),& + sinph1ph2(l,k),sinph1ph2(k,l) + enddo + enddo + write(iout,*) "ethetai",ethetai + endif + do m=1,ntheterm2_nucl + do k=1,nsingle_nucl + aux=bbthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)& + +ccthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)& + +ddthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)& + +eethet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k) + ethetai=ethetai+sinkt(m)*aux + dethetai=dethetai+0.5d0*m*aux*coskt(m) + dephii=dephii+k*sinkt(m)*(& + ccthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)-& + bbthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)) + dephii1=dephii1+k*sinkt(m)*(& + eethet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)-& + ddthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k)) + if (lprn) & + write (iout,*) "m",m," k",k," bbthet",& + bbthet_nucl(k,m,ityp1,ityp2,ityp3)," ccthet",& + ccthet_nucl(k,m,ityp1,ityp2,ityp3)," ddthet",& + ddthet_nucl(k,m,ityp1,ityp2,ityp3)," eethet",& + eethet_nucl(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai + enddo + enddo + if (lprn) & + write(iout,*) "ethetai",ethetai + do m=1,ntheterm3_nucl + do k=2,ndouble_nucl + do l=1,k-1 + aux=ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+& + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l) + ethetai=ethetai+sinkt(m)*aux + dethetai=dethetai+0.5d0*m*coskt(m)*aux + dephii=dephii+l*sinkt(m)*(& + -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-& + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + dephii1=dephii1+(k-l)*sinkt(m)*( & + -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+& + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + if (lprn) then + write (iout,*) "m",m," k",k," l",l," ffthet", & + ffthet_nucl(l,k,m,ityp1,ityp2,ityp3), & + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ggthet",& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3),& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai + write (iout,*) cosph1ph2(l,k)*sinkt(m), & + cosph1ph2(k,l)*sinkt(m),& + sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m) + endif + enddo + enddo + enddo +10 continue + if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & + i,theta(i)*rad2deg,phii*rad2deg, & + phii1*rad2deg,ethetai + etheta_nucl=etheta_nucl+ethetai + print *,i,"partial sum",etheta_nucl + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii + if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1 + gloc(nphi+i-2,icg)=wang_nucl*dethetai + enddo + return + end subroutine ebend_nucl + !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module energy