X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=2c663fb729e074ddb5016e93237697e0a062d686;hb=c8b71b34ec6b864cc5db11cc5c0d95349a9b5cef;hp=94ea0332c47dc90b4e021d54e01baff9665dd1b1;hpb=5bb791af1e8990b66bc0543222e13e789a3f69aa;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 94ea033..2c663fb 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 /.false./ 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 @@ -3122,8 +3171,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 +3197,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 +3212,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 +3281,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 +3424,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 +3439,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 +3453,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 +3475,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 +3593,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 +3985,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 +4579,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 +5288,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 +5380,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 +5535,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 +5544,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 @@ -6435,7 +6498,7 @@ c---------------------------------------------------------------------------- include 'COMMON.GEO' logical swap double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2), - & auxvec1(2),auxvec2(1),auxmat1(2,2) + & auxvec1(2),auxvec2(2),auxmat1(2,2) logical lprn common /kutas/ lprn CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC @@ -6640,14 +6703,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 +6825,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