merge...

Merge branch 'prerelease-3.2.1' of mmka:unres into prerelease-3.2.1

Conflicts:

	source/wham/src/energy_p_new.F

diff --combined source/wham/src/energy_p_new.F
index 12f60f5,b6540fb..5b36ced
--- 1/source/wham/src/energy_p_new.F
--- 2/source/wham/src/energy_p_new.F
+++ b/source/wham/src/energy_p_new.F
@@@ -107,7 -107,7 +107,7 @@@ c      write (iout,*) "ft(6)",fact(6),
        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 +116,7 @@@
        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
@@@ -154,6 -154,7 +154,7 @@@
        energia(19)=esccor
        energia(20)=edihcnstr
        energia(21)=evdw_t
+ c      if (dyn_ss) call dyn_set_nss
  c detecting NaNQ
  #ifdef ISNAN
  #ifdef AIX
@@@ -770,6 -771,7 +771,7 @@@
        include 'COMMON.ENEPS'
        include 'COMMON.IOUNITS'
        include 'COMMON.CALC'
+       include 'COMMON.SBRIDGE'
        logical lprn
        common /srutu/icall
        integer icant
@@@ -800,6 -802,21 +802,21 @@@ C Calculate SC interaction energy
  C
          do iint=1,nint_gr(i)
            do j=istart(i,iint),iend(i,iint)
+ C in case of diagnostics    write (iout,*) "TU SZUKAJ",i,j,dyn_ss_mask(i),dyn_ss_mask(j)
+ C /06/28/2013 Adasko: In case of dyn_ss - dynamic disulfide bond
+ C formation no electrostatic interactions should be calculated. If it
+ C would be allowed NaN would appear
+             IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+ C /06/28/2013 Adasko: dyn_ss_mask is logical statement wheather this Cys
+ C residue can or cannot form disulfide bond. There is still bug allowing
+ C Cys...Cys...Cys bond formation
+               call dyn_ssbond_ene(i,j,evdwij)
+ C /06/28/2013 Adasko: dyn_ssbond_ene is dynamic SS bond foration energy
+ C function in ssMD.F
+               evdw=evdw+evdwij
+ c              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+ c     &                        'evdw',i,j,evdwij,' ss'
+             ELSE
              ind=ind+1
              itypj=itype(j)
              dscj_inv=vbld_inv(j+nres)
@@@ -866,6 -883,7 +883,7 @@@ c--------------------------------------
  c            write (iout,*) "i",i," j",j," itypi",itypi," itypj",itypj,
  c     &         " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
  c     &         aux*e2/eps(itypi,itypj)
+ c       write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
              if (lprn) then
              sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
              epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@@ -889,6 -907,7 +907,7 @@@ C Calculate the radial part of the grad
  C Calculate angular part of the gradient.
              call sc_grad
              endif
+             ENDIF    ! dyn_ss
            enddo      ! j
          enddo        ! iint
        enddo          ! i
@@@ -2897,9 -2916,12 +2916,12 @@@ c        write (iout,*) "i",i," ii",ii,
  c     &    dhpb(i),dhpb1(i),forcon(i)
  C 24/11/03 AL: SS bridges handled separately because of introducing a specific
  C    distance and angle dependent SS bond potential.
+         if (.not.dyn_ss .and. i.le.nss) then
+ C 15/02/13 CC dynamic SSbond - additional check
          if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
            call ssbond_ene(iii,jjj,eij)
            ehpb=ehpb+2*eij
+          endif
  cd          write (iout,*) "eij",eij
          else if (ii.gt.nres .and. jj.gt.nres) then
  c Restraints from contact prediction
@@@ -3032,11 -3054,12 +3054,12 @@@
        deltat12=om2-om1+2.0d0
        cosphi=om12-om1*om2
        eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
-      &  +akct*deltad*deltat12
+      &  +akct*deltad*deltat12+ebr
+ c     &  +akct*deltad*deltat12
       &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
- c      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
- c     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
- c     &  " deltat12",deltat12," eij",eij 
+       write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
+      &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
+      &  " deltat12",deltat12," eij",eij,"ebr",ebr
        ed=2*akcm*deltad+akct*deltat12
        pom1=akct*deltad
        pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
@@@ -3081,6 -3104,7 +3104,7 @@@
        include 'COMMON.FFIELD'
        include 'COMMON.CONTROL'
        double precision u(3),ud(3)
+       logical :: lprn=.false.
        estr=0.0d0
        do i=nnt+1,nct
          diff = vbld(i)-vbldp0
@@@ -3100,8 -3124,9 +3124,14 @@@
            nbi=nbondterm(iti)
            if (nbi.eq.1) then
              diff=vbld(i+nres)-vbldsc0(1,iti)
++<<<<<<< HEAD
 +c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
 +c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
++=======
+             if (lprn)
+      &      write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+      &      AKSC(1,iti),AKSC(1,iti)*diff*diff
++>>>>>>> aee20d3590dc2913e3a9a4308ce5da7787993a66
              estr=estr+0.5d0*AKSC(1,iti)*diff*diff
              do j=1,3
                gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@@ -3130,8 -3155,9 +3160,9 @@@
                usum=usum+uprod1
                usumsqder=usumsqder+ud(j)*uprod2
              enddo
- c            write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
- c     &      AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
+             if (lprn)
+      &      write (iout,*) i,iti,vbld(i+nres),(vbldsc0(j,iti),
+      &      AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi)
              estr=estr+uprod/usum
              do j=1,3
               gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres)
@@@ -4003,7 -4029,8 +4034,8 @@@ c        sumene = enesc(x,xx,yy,zz,cost
          escloc = escloc + sumene
  c        write (2,*) "escloc",escloc
          if (.not. calc_grad) goto 1
- #ifdef DEBUG
+ 
+ #ifdef DEBUG2
  C
  C This section to check the numerical derivatives of the energy of ith side
  C chain in xx, yy, zz, and theta. Use the -DDEBUG compiler option or insert

diff --combined source/wham/src/wham_calc1.F
index 57f3fcf,676d4f4..61b9504
--- 1/source/wham/src/wham_calc1.F
--- 2/source/wham/src/wham_calc1.F
+++ b/source/wham/src/wham_calc1.F
@@@ -94,6 -94,7 +94,6 @@@ c      parameter (MaxHdim=200000
        character*128 nazwa
        integer ilen
        external ilen
 - 
        write(licz2,'(bz,i2.2)') islice
        nbin1 = 1.0d0/delta
        write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
@@@ -315,7 -316,7 +315,7 @@@ c              write (iout,*) 1.0d0/(be
              etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
       &      +wvdwpp*evdw1
       &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
       &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
       &      +ft(2)*wturn3*eello_turn3
       &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
@@@ -325,7 -326,7 +325,7 @@@
              etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
       &      +ft(1)*welec*(ees+evdw1)
       &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
       &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
       &      +ft(2)*wturn3*eello_turn3
       &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
@@@ -674,7 -675,7 +674,7 @@@ c              write (iout,*) 1.0d0/(be
              etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
       &      +wvdwpp*evdw1
       &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
       &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
       &      +ft(2)*wturn3*eello_turn3
       &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
@@@ -684,7 -685,7 +684,7 @@@
              etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
       &      +ft(1)*welec*(ees+evdw1)
       &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
       &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
       &      +ft(2)*wturn3*eello_turn3
       &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
@@@ -1012,7 -1013,7 +1012,7 @@@ c            write (iout,*) ib," PotEmi
              etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
       &      +wvdwpp*evdw1
       &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
       &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
       &      +ft(2)*wturn3*eello_turn3
       &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
@@@ -1035,7 -1036,7 +1035,7 @@@
              etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
       &      +ft(1)*welec*(ees+evdw1)
       &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
       &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
       &      +ft(2)*wturn3*eello_turn3
       &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr