X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=a65c91cff30a024ecc99b1f49a77dd4807c2d486;hb=03b8a72f0445bc2539b8d3c286fb15ca24c0ff16;hp=c588bee9723d874787dd199c98bbce2e093bedeb;hpb=462b4b0ab826c82da8fbe9e4425f5329306b1952;p=unres.git diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index c588bee..a65c91c 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -107,7 +107,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees & +wvdwpp*evdw1 & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d @@ -116,7 +116,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) & +wang*ebe+wtor*fact(1)*etors+wscloc*escloc - & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 + & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5 & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d @@ -124,7 +124,6 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t #endif energia(0)=etot energia(1)=evdw -c call enerprint(energia(0),frac) #ifdef SCP14 energia(2)=evdw2-evdw2_14 energia(17)=evdw2_14 @@ -229,9 +228,11 @@ C & +wturn3*fact(2)*gel_loc_turn3(i) & +wturn6*fact(5)*gel_loc_turn6(i) & +wel_loc*fact(2)*gel_loc_loc(i) - & +wsccor*fact(1)*gsccor_loc(i) +c & +wsccor*fact(1)*gsccor_loc(i) +c ROZNICA Z WHAMem enddo endif + if (dyn_ss) call dyn_set_nss return end C------------------------------------------------------------------------ @@ -360,6 +361,14 @@ C integer icant external icant cd print *,'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon +c ROZNICA DODANE Z WHAM +c do i=1,210 +c do j=1,2 +c eneps_temp(j,i)=0.0d0 +c enddo +c enddo +cROZNICA + evdw=0.0D0 evdw_t=0.0d0 do i=iatsc_s,iatsc_e @@ -393,6 +402,11 @@ c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj e2=fac*bb(itypi,itypj) evdwij=e1+e2 ij=icant(itypi,itypj) +c ROZNICA z WHAM +c eneps_temp(1,ij)=eneps_temp(1,ij)+e1/dabs(eps0ij) +c eneps_temp(2,ij)=eneps_temp(2,ij)+e2/eps0ij +c + cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0) cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj) cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)') @@ -750,10 +764,12 @@ C include 'COMMON.INTERACT' include 'COMMON.IOUNITS' include 'COMMON.CALC' + include 'COMMON.SBRIDGE' logical lprn common /srutu/icall integer icant external icant + logical energy_dec /.true./ c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon evdw=0.0D0 evdw_t=0.0d0 @@ -776,6 +792,38 @@ C Calculate SC interaction energy. C do iint=1,nint_gr(i) do j=istart(i,iint),iend(i,iint) + IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN + +c write(iout,*) "PRZED ZWYKLE", evdwij + call dyn_ssbond_ene(i,j,evdwij) +c write(iout,*) "PO ZWYKLE", evdwij + + evdw=evdw+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) +C search over all next residues + 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) +c if(evdwij_przed_tri.ne.evdwij) then +c write (iout,*) "TRI:", evdwij, evdwij_przed_tri +c endif + +c write(iout,*) "PO TRI", evdwij +C call the energy function that removes the artifical triple disulfide +C bond the soubroutine is located in ssMD.F + evdw=evdw+evdwij + if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') + & 'evdw',i,j,evdwij,'tss' + endif!dyn_ss_mask(k) + enddo! k + ELSE ind=ind+1 itypj=iabs(itype(j)) if (itypj.eq.ntyp1) cycle @@ -864,6 +912,7 @@ C Calculate the radial part of the gradient C Calculate angular part of the gradient. call sc_grad endif + ENDIF ! dyn_ss enddo ! j enddo ! iint enddo ! i @@ -2880,11 +2929,41 @@ C iii and jjj point to the residues for which the distance is assigned. endif C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C & iabs(itype(jjj)).eq.1) then +C call ssbond_ene(iii,jjj,eij) +C ehpb=ehpb+2*eij +C else + if (.not.dyn_ss .and. i.le.nss) then + if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. + & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij - else + endif !ii.gt.neres + else if (ii.gt.nres .and. jj.gt.nres) then +c Restraints from contact prediction + dd=dist(ii,jj) + if (constr_dist.eq.11) then +C ehpb=ehpb+fordepth(i)**4.0d0 +C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd +C write(iout,*) ehpb,"atu?" +C ehpb,"tu?" +C fac=fordepth(i)**4.0d0 +C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd + else !constr_dist.eq.11 + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "beta nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else !dhpb(i).gt.0.00 + C Calculate the distance between the two points and its difference from the C target distance. dd=dist(ii,jj) @@ -2897,6 +2976,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd + endif !dhpb(i).gt.0 + endif cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd, cd & ' waga=',waga,' fac=',fac do j=1,3 @@ -2911,6 +2992,53 @@ C Cartesian gradient in the SC vectors (ghpbx). ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) enddo endif + else !ii.gt.nres +C write(iout,*) "before" + dd=dist(ii,jj) +C write(iout,*) "after",dd + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) +C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd +C print *,ehpb,"tu?" +C write(iout,*) ehpb,"btu?", +C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else + rdis=dd-dhpb(i) +C Get the force constant corresponding to this distance. + waga=forcon(i) +C Calculate the contribution to energy. + ehpb=ehpb+waga*rdis*rdis +c write (iout,*) "alpha reg",dd,waga*rdis*rdis +C +C Evaluate gradient. +C + fac=waga*rdis/dd + endif + endif + do j=1,3 + ggg(j)=fac*(c(j,jj)-c(j,ii)) + enddo +cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3) +C If this is a SC-SC distance, we need to calculate the contributions to the +C Cartesian gradient in the SC vectors (ghpbx). + if (iii.lt.ii) then + do j=1,3 + ghpbx(j,iii)=ghpbx(j,iii)-ggg(j) + ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j) + enddo + endif do j=iii,jjj-1 do k=1,3 ghpbc(k,j)=ghpbc(k,j)+ggg(k) @@ -2918,7 +3046,7 @@ C Cartesian gradient in the SC vectors (ghpbx). enddo endif enddo - ehpb=0.5D0*ehpb + if (constr_dist.ne.11) ehpb=0.5D0*ehpb return end C-------------------------------------------------------------------------- @@ -3122,8 +3250,8 @@ C & delthe0,sig0inv,sigtc,sigsqtc,delthec,it double precision y(2),z(2) delta=0.02d0*pi - time11=dexp(-2*time) - time12=1.0d0 +c time11=dexp(-2*time) +c time12=1.0d0 etheta=0.0D0 c write (iout,*) "nres",nres c write (*,'(a,i2)') 'EBEND ICG=',icg @@ -3148,8 +3276,8 @@ C Zero the energy function and its derivative at 0 or pi. if (i.gt.3 .and. itype(i-2).ne.ntyp1) then #ifdef OSF phii=phi(i) - icrc=0 - call proc_proc(phii,icrc) +c icrc=0 +c call proc_proc(phii,icrc) if (icrc.eq.1) phii=150.0 #else phii=phi(i) @@ -3163,8 +3291,8 @@ C Zero the energy function and its derivative at 0 or pi. if (i.lt.nres .and. itype(i).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) - icrc=0 - call proc_proc(phii1,icrc) +c icrc=0 +c call proc_proc(phii1,icrc) if (icrc.eq.1) phii1=150.0 phii1=pinorm(phii1) z(1)=cos(phii1) @@ -3232,7 +3360,7 @@ c & rad2deg*phii,rad2deg*phii1,ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) - 1215 continue +c 1215 continue enddo C Ufff.... We've done all this!!! return @@ -3375,7 +3503,9 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end - if (itype(i-1).eq.ntyp1) cycle +c if (itype(i-1).eq.ntyp1) cycle + if ((itype(i-1).eq.ntyp1).or.(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 dethetai=0.0d0 @@ -3388,7 +3518,7 @@ CC Ta zmina jest niewlasciwa coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.ntyp1) then + if (i.gt.3 .and. itype(i-3).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -3402,13 +3532,14 @@ CC Ta zmina jest niewlasciwa enddo else phii=0.0d0 - ityp1=nthetyp+1 +c ityp1=nthetyp+1 do k=1,nsingle + ityp1=ithetyp((itype(i-2))) cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo 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 @@ -3423,7 +3554,8 @@ CC Ta zmina jest niewlasciwa enddo else phii1=0.0d0 - ityp3=nthetyp+1 +c ityp3=nthetyp+1 + ityp3=ithetyp((itype(i))) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 @@ -3540,7 +3672,8 @@ c call flush(iout) etheta=etheta+ethetai if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 - gloc(nphi+i-2,icg)=wang*dethetai +c gloc(nphi+i-2,icg)=wang*dethetai + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo return end @@ -3931,7 +4064,8 @@ c write (2,*) "xx",xx," yy",yy," zz",zz Cc diagnostics - remove later xx1 = dcos(alph(2)) yy1 = dsin(alph(2))*dcos(omeg(2)) - zz1 = -dsin(alph(2))*dsin(omeg(2)) +c zz1 = -dsin(alph(2))*dsin(omeg(2)) + zz1 = -dsign(1.0d0,itype(i))*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 @@ -4524,7 +4658,7 @@ c lprn=.true. c write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor 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).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle esccor_ii=0.0D0 isccori=isccortyp(itype(i-2)) isccori1=isccortyp(itype(i-1)) @@ -5233,7 +5367,11 @@ C--------------------------------------------------------------------------- & auxmat(2,2) iti1 = itortyp(itype(i+1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1)) + if (itype(j).le.ntyp) then + itj1 = itortyp(itype(j+1)) + else + itj1=ntortyp+1 + endif else itj1=ntortyp+1 endif @@ -5321,14 +5459,16 @@ cd if (i.ne.2 .or. j.ne.4 .or. k.ne.3 .or. l.ne.5) return enddo if (l.eq.j+1) then C parallel orientation of the two CA-CA-CA frames. - if (i.gt.1) then +c if (i.gt.1) then + if (i.gt.1 .and. itype(i).le.ntyp) then iti=itortyp(itype(i)) else iti=ntortyp+1 endif itk1=itortyp(itype(k+1)) itj=itortyp(itype(j)) - if (l.lt.nres-1) then +c if (l.lt.nres-1) then + if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then itl1=itortyp(itype(l+1)) else itl1=ntortyp+1 @@ -5474,7 +5614,8 @@ C Calculate the Cartesian derivatives of the vectors. C End vectors else C Antiparallel orientation of the two CA-CA-CA frames. - if (i.gt.1) then +c if (i.gt.1) then + if (i.gt.1 .and. itype(i).le.ntyp) then iti=itortyp(itype(i)) else iti=ntortyp+1 @@ -5482,7 +5623,8 @@ C Antiparallel orientation of the two CA-CA-CA frames. itk1=itortyp(itype(k+1)) itl=itortyp(itype(l)) itj=itortyp(itype(j)) - if (j.lt.nres-1) then +c if (j.lt.nres-1) then + if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then itj1=itortyp(itype(j+1)) else itj1=ntortyp+1 @@ -6640,14 +6782,16 @@ C C 4/7/01 AL Component s1 was removed, because it pertains to the respective C energy moment and not to the cluster cumulant. iti=itortyp(itype(i)) - if (j.lt.nres-1) then +c if (j.lt.nres-1) then + if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then itj1=itortyp(itype(j+1)) else itj1=ntortyp+1 endif itk=itortyp(itype(k)) itk1=itortyp(itype(k+1)) - if (l.lt.nres-1) then +c if (l.lt.nres-1) then + if (l.lt.nres-1 .and. itype(l+1).le.ntyp) then itl1=itortyp(itype(l+1)) else itl1=ntortyp+1 @@ -6760,13 +6904,15 @@ C energy moment and not to the cluster cumulant. cd write (2,*) 'eello_graph4: wturn6',wturn6 iti=itortyp(itype(i)) itj=itortyp(itype(j)) - if (j.lt.nres-1) then +c if (j.lt.nres-1) then + if (j.lt.nres-1 .and. itype(j+1).le.ntyp) then itj1=itortyp(itype(j+1)) else itj1=ntortyp+1 endif itk=itortyp(itype(k)) - if (k.lt.nres-1) then +c if (k.lt.nres-1) then + if (k.lt.nres-1 .and. itype(k+1).le.ntyp) then itk1=itortyp(itype(k+1)) else itk1=ntortyp+1