Merge branch 'devel' of mmka:unres into czarek

Conflicts:
	bin/unres/MD/unres_ifort_MPICH_GAB.exe
	bin/unres/MD/unres_ifort_single_GAB.exe
	source/unres/src_MD/parmread.F

diff --combined bin/unres/MD/unres_ifort_MPICH_GAB.exe
index 886e8d6,8099588..0000000
deleted file mode 100755,100755
Binary files differ

diff --combined bin/unres/MD/unres_ifort_single_GAB.exe
index d08cbf4,7107570..0000000
deleted file mode 100755,100755
Binary files differ

diff --combined source/unres/src_MD/CMakeLists.txt
index ceac982,b041904..94d735a
--- 1/source/unres/src_MD/CMakeLists.txt
--- 2/source/unres/src_MD/CMakeLists.txt
+++ b/source/unres/src_MD/CMakeLists.txt
@@@ -78,11 -78,12 +78,13 @@@ set(UNRES_MD_SRC
        timing.F
        thread.F 
        unres.F
 +      ssMD.F
  )
  
  if(Fortran_COMPILER_NAME STREQUAL "ifort")
    set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng.f ) 
+ elseif(Fortran_COMPILER_NAME STREQUAL "mpif90")
+   set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng.f )
  else()
    set(UNRES_MD_SRC0 ${UNRES_MD_SRC0} prng_32.F )
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
@@@ -150,11 -151,11 +152,11 @@@ if (Fortran_COMPILER_NAME STREQUAL "ifo
    #set(FFLAGS3 "-c -w -O3 -ipo -ipo_obj -opt_report" )
    set(FFLAGS3 "-w -ipo " )
  elseif (Fortran_COMPILER_NAME STREQUAL "gfortran")
-   set(FFLAGS0 "-I. " ) 
-   set(FFLAGS1 "-g -I. " ) 
-   set(FFLAGS2 "-I. ")
+   set(FFLAGS0 "-std=legacy -I. " ) 
+   set(FFLAGS1 "-std=legacy -g -I. " ) 
+   set(FFLAGS2 "-std=legacy -I. ")
    #set(FFLAGS3 "-c -w -O3 -ipo -ipo_obj -opt_report" )
-   set(FFLAGS3 "-I. " )
+   set(FFLAGS3 "-std=legacy -I. " )
  endif (Fortran_COMPILER_NAME STREQUAL "ifort")
  
  
@@@ -261,9 -262,6 +263,6 @@@ FILE(APPEND ${CINFO
         return 
         end ")
  
- #FILE(APPEND ${CMAKE_CURRENT_BINARY_DIR}/cinfo.f
- #     CINFO_FORMAT(CPPFLAGS)
- #)
  # add include path
  set_property(SOURCE ${CMAKE_CURRENT_BINARY_DIR}/cinfo.f PROPERTY COMPILE_FLAGS "${FFLAGS0} -I${CMAKE_CURRENT_SOURCE_DIR}")
  

diff --combined source/unres/src_MD/Makefile_single_ifort
index a9fdabf,245206b..c651e39
--- 1/source/unres/src_MD/Makefile_single_ifort
--- 2/source/unres/src_MD/Makefile_single_ifort
+++ b/source/unres/src_MD/Makefile_single_ifort
@@@ -37,13 -37,13 +37,13 @@@ object = unres.o arcos.o cartprint.o ch
          cored.o rmdd.o geomout.o readpdb.o regularize.o thread.o fitsq.o mcm.o \
          mc.o bond_move.o refsys.o check_sc_distr.o check_bond.o contact.o djacob.o \
          eigen.o blas.o add.o entmcm.o minim_mcmf.o \
-         MP.o compare_s1.o prng_32.o \
+         MP.o compare_s1.o prng.o \
          banach.o rmsd.o elecont.o dihed_cons.o \
          sc_move.o local_move.o \
          intcartderiv.o lagrangian_lesyng.o\
        stochfric.o kinetic_lesyng.o MD_A-MTS.o moments.o int_to_cart.o \
          surfatom.o sort.o muca_md.o MREMD.o rattle.o gauss.o energy_split-sep.o \
 -        q_measure.o gnmr1.o test.o
 +        q_measure.o gnmr1.o test.o ssMD.o
  
  GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN \
        -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC

diff --combined source/unres/src_MD/energy_p_new_barrier.F
index 2fe5fdb,981a461..3001a22
--- 1/source/unres/src_MD/energy_p_new_barrier.F
--- 2/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@@ -131,11 -131,6 +131,11 @@@
  C Calculate electrostatic (H-bonding) energy of the main chain.
  C
    107 continue
 +cmc
 +cmc Sep-06: egb takes care of dynamic ss bonds too
 +cmc
 +      if (dyn_ss) call dyn_set_nss
 +
  c      print *,"Processor",myrank," computed USCSC"
  #ifdef TIMING
  #ifdef MPI
@@@ -428,14 -423,14 +428,14 @@@ cMS$ATTRIBUTES C ::  proc_pro
  #ifdef SPLITELE
        etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
       & +wang*ebe+wtor*etors+wscloc*escloc
-      & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
+      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
       & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
       & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
       & +wbond*estr+Uconst+wsccor*esccor
  #else
        etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
       & +wang*ebe+wtor*etors+wscloc*escloc
-      & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
+      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
       & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
       & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
       & +wbond*estr+Uconst+wsccor*esccor
@@@ -1057,25 -1052,25 +1057,25 @@@ C--------------------------------------
       &  edihcnstr,ebr*nss,
       &  Uconst,etot
     10 format (/'Virtual-chain energies:'//
-      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
-      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
-      & 'EES=   ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/
-      & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/
-      & 'ESTR=  ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/
-      & 'EBE=   ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/
-      & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
-      & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
-      & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
-      & 'EHPB=  ',1pE16.6,' WEIGHT=',1pD16.6,
+      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/
+      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/
+      & 'EES=   ',1pE16.6,' WEIGHT=',1pE16.6,' (p-p)'/
+      & 'EVDWPP=',1pE16.6,' WEIGHT=',1pE16.6,' (p-p VDW)'/
+      & 'ESTR=  ',1pE16.6,' WEIGHT=',1pE16.6,' (stretching)'/
+      & 'EBE=   ',1pE16.6,' WEIGHT=',1pE16.6,' (bending)'/
+      & 'ESC=   ',1pE16.6,' WEIGHT=',1pE16.6,' (SC local)'/
+      & 'ETORS= ',1pE16.6,' WEIGHT=',1pE16.6,' (torsional)'/
+      & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
+      & 'EHPB=  ',1pE16.6,' WEIGHT=',1pE16.6,
       & ' (SS bridges & dist. cnstr.)'/
-      & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-      & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-      & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-      & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/
-      & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/
-      & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/
-      & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
-      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
+      & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+      & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+      & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+      & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
+      & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
+      & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+      & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+      & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
       & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
       & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
       & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
@@@ -1557,7 -1552,6 +1557,7 @@@
        include 'COMMON.IOUNITS'
        include 'COMMON.CALC'
        include 'COMMON.CONTROL'
 +      include 'COMMON.SBRIDGE'
        logical lprn
        evdw=0.0D0
  ccccc      energy_dec=.false.
@@@ -1586,10 -1580,6 +1586,10 @@@ 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
 +              call dyn_ssbond_ene(i,j,evdwij)
 +              evdw=evdw+evdwij
 +            ELSE
              ind=ind+1
              itypj=itype(j)
  c            dscj_inv=dsc_inv(itypj)
@@@ -1699,7 -1689,6 +1699,7 @@@ C Calculate angular part of the gradien
  #else
              call sc_grad
  #endif
 +            ENDIF    ! dyn_ss            
            enddo      ! j
          enddo        ! iint
        enddo          ! i
@@@ -4257,8 -4246,8 +4257,8 @@@
        include 'COMMON.IOUNITS'
        dimension ggg(3)
        ehpb=0.0D0
 -cd      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
 -cd      write(iout,*)'link_start=',link_start,' link_end=',link_end
 +      write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
 +      write(iout,*)'link_start=',link_start,' link_end=',link_end
        if (link_end.eq.0) return
        do i=link_start,link_end
  C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
@@@ -4280,7 -4269,7 +4280,7 @@@ C    distance and angle dependent SS bo
          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
 -cd          write (iout,*) "eij",eij
 +          write (iout,*) "eij",eij
          else if (ii.gt.nres .and. jj.gt.nres) then
  c Restraints from contact prediction
            dd=dist(ii,jj)
@@@ -4418,9 -4407,8 +4418,9 @@@ c      dscj_inv=dsc_inv(itypj
        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
       &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 +     &  +ss_depth
  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 

diff --combined source/unres/src_MD/initialize_p.F
index 01628ec,ec978e7..5585cae
--- 1/source/unres/src_MD/initialize_p.F
--- 2/source/unres/src_MD/initialize_p.F
+++ b/source/unres/src_MD/initialize_p.F
@@@ -345,7 -345,6 +345,7 @@@ cd    write (iout,*) 'ns=',ns,' nss=',n
  cd   &   (ihpb(i),jhpb(i),i=1,nss)
        do i=nnt,nct-1
          scheck=.false.
 +        if (dyn_ss) goto 10
          do ii=1,nss
            if (ihpb(ii).eq.i+nres) then
              scheck=.true.
@@@ -1376,6 -1375,7 +1376,7 @@@ c--------------------------------------
        include 'COMMON.IOUNITS'
        include 'COMMON.SETUP'
        include 'COMMON.CONTROL'
+ c      write(2,*)"hpb_partition: nhpb=",nhpb
  #ifdef MPI
        call int_bounds(nhpb,link_start,link_end)
        if (.not. out1file) 
@@@ -1387,5 -1387,6 +1388,6 @@@
        link_start=1
        link_end=nhpb
  #endif
+ c      write(2,*)"hpb_partition: link_start=",nhpb," link_end=",link_end
        return
        end

diff --combined source/unres/src_MD/readrtns.F
index 0d8fee2,5b74f87..a11f4fb
--- 1/source/unres/src_MD/readrtns.F
--- 2/source/unres/src_MD/readrtns.F
+++ b/source/unres/src_MD/readrtns.F
@@@ -860,36 -860,12 +860,36 @@@ C 12/1/95 Added weight for the multi-bo
        call reada(weightcard,"V2SS",v2ss,7.61d0)
        call reada(weightcard,"V3SS",v3ss,13.7d0)
        call reada(weightcard,"EBR",ebr,-5.50D0)
 +      dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 +      do i=1,maxres
 +        dyn_ss_mask(i)=.false.
 +      enddo
 +      do i=1,maxres-1
 +        do j=i+1,maxres
 +          dyn_ssbond_ij(i,j)=1.0d300
 +        enddo
 +      enddo
 +      call reada(weightcard,"HT",Ht,0.0D0)
 +      if (dyn_ss) then
 +        ss_depth=ebr/wsc-0.25*eps(1,1)
 +        Ht=Ht/wsc-0.25*eps(1,1)
 +        akcm=akcm*wstrain/wsc
 +        akth=akth*wstrain/wsc
 +        akct=akct*wstrain/wsc
 +        v1ss=v1ss*wstrain/wsc
 +        v2ss=v2ss*wstrain/wsc
 +        v3ss=v3ss*wstrain/wsc
 +      else
 +        ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 +      endif
 +
        if(me.eq.king.or..not.out1file) then
         write (iout,*) "Parameters of the SS-bond potential:"
         write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
       & " AKCT",akct
         write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 -       write (iout,*) "EBR",ebr
 +       write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 +       write (iout,*)" HT",Ht
         print *,'indpdb=',indpdb,' pdbref=',pdbref
        endif
        if (indpdb.gt.0 .or. pdbref) then
@@@ -917,6 -893,9 +917,9 @@@ c        print *,'Finished reading pdb 
          call contact(.false.,ncont_ref,icont_ref,co)
  
          if (sideadd) then 
+ C Following 2 lines for diagnostics; comment out if not needed
+          write (iout,*) "Before sideadd"
+          call intout
           if(me.eq.king.or..not.out1file)
       &    write(iout,*)'Adding sidechains'
           maxsi=1000
@@@ -933,7 -912,12 +936,12 @@@
       &              i,' after ',nsi,' trials'
            endif
           enddo
+ C 10/03/12 Adam: Recalculate coordinates with new side chain positions
+          call chainbuild
          endif  
+ C Following 2 lines for diagnostics; comment out if not needed
+         write (iout,*) "After sideadd"
+         call intout
        endif
        if (indpdb.eq.0) then
  C Read sequence if not taken from the pdb file.
@@@ -1241,21 -1225,6 +1249,21 @@@ C Generate distance constraints, if th
          write (iout,'(/a,i3,a)') 
       &  'The chain contains',ns,' disulfide-bridging cysteines.'
          write (iout,'(20i4)') (iss(i),i=1,ns)
 +       if (dyn_ss) then
 +          write(iout,*)"Running with dynamic disulfide-bond formation"
 +          do i=nss+1,nhpb
 +            ihpb(i-nss)=ihpb(i)
 +            jhpb(i-nss)=jhpb(i)
 +            forcon(i-nss)=forcon(i)
 +            dhpb(i-nss)=dhpb(i)
 +          enddo
 +          nhpb=nhpb-nss
 +          nss=0
 +          call hpb_partition
 +          do i=1,ns
 +            dyn_ss_mask(iss(i))=.true.
 +          enddo
 +       else
          write (iout,'(/a/)') 'Pre-formed links are:' 
        do i=1,nss
          i1=ihpb(i)-nres
@@@ -1268,7 -1237,6 +1276,7 @@@
       &    ebr,forcon(i)
        enddo
        write (iout,'(a)')
 +       endif
        endif
        if (i2ndstr.gt.0) call secstrp2dihc
  c      call geom_to_var(nvar,x)