dynamic SS from old Maurizio's code
[unres.git] / source / unres / src_MD / readrtns.F
index 4d8eb72..04b809c 100644 (file)
@@ -138,7 +138,7 @@ C Set up the time limit (caution! The time must be input in minutes!)
       call readi(controlcard,'MAXGEN',maxgen,10000)
       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
       call readi(controlcard,"KDIAG",kdiag,0)
-      call readi(controlcard,"RESCALE_MODE",rescale_mode,1)
+      call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
      & write (iout,*) "RESCALE_MODE",rescale_mode
       split_ene=index(controlcard,'SPLIT_ENE').gt.0
@@ -860,12 +860,36 @@ C 12/1/95 Added weight for the multi-body term WCORR
       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
@@ -1217,6 +1241,21 @@ C Generate distance constraints, if the PDB structure is to be regularized.
         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
@@ -1229,6 +1268,7 @@ C Generate distance constraints, if the PDB structure is to be regularized.
      &    ebr,forcon(i)
        enddo
        write (iout,'(a)')
+       endif
       endif
       if (i2ndstr.gt.0) call secstrp2dihc
 c      call geom_to_var(nvar,x)
@@ -1813,6 +1853,9 @@ c----------------------------------------------------------------------------
       call readi(minimcard,'MINFUN',minfun,maxmin)
       call reada(minimcard,'TOLF',tolf,1.0D-2)
       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
+      print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
+      print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
+      print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
       write (iout,'(/80(1h*)/20x,a/80(1h*))') 
      &         'Options in energy minimization:'
       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
@@ -2083,38 +2126,38 @@ C Get parameter filenames and open the parameter files.
       open (isidep,file=sidename,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
-     &  readonly)
+     &action='read')
        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
 C Get parameter filenames and open the parameter files.
       call getenv_loc('BONDPAR',bondname)
-      open (ibond,file=bondname,status='old',readonly)
+      open (ibond,file=bondname,status='old',action='read')
       call getenv_loc('THETPAR',thetname)
-      open (ithep,file=thetname,status='old',readonly)
+      open (ithep,file=thetname,status='old',action='read')
 #ifndef CRYST_THETA
       call getenv_loc('THETPARPDB',thetname_pdb)
       print *,"thetname_pdb ",thetname_pdb
-      open (ithep_pdb,file=thetname_pdb,status='old',readonly)
+      open (ithep_pdb,file=thetname_pdb,status='old',action='read')
       print *,ithep_pdb," opened"
 #endif
       call getenv_loc('ROTPAR',rotname)
-      open (irotam,file=rotname,status='old',readonly)
+      open (irotam,file=rotname,status='old',action='read')
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
-      open (irotam_pdb,file=rotname_pdb,status='old',readonly)
+      open (irotam_pdb,file=rotname_pdb,status='old',action='read')
 #endif
       call getenv_loc('TORPAR',torname)
-      open (itorp,file=torname,status='old',readonly)
+      open (itorp,file=torname,status='old',action='read')
       call getenv_loc('TORDPAR',tordname)
-      open (itordp,file=tordname,status='old',readonly)
+      open (itordp,file=tordname,status='old',action='read')
       call getenv_loc('SCCORPAR',sccorname)
-      open (isccor,file=sccorname,status='old',readonly)
+      open (isccor,file=sccorname,status='old',action='read')
       call getenv_loc('FOURIER',fouriername)
-      open (ifourier,file=fouriername,status='old',readonly)
+      open (ifourier,file=fouriername,status='old',action='read')
       call getenv_loc('ELEPAR',elename)
-      open (ielep,file=elename,status='old',readonly)
+      open (ielep,file=elename,status='old',action='read')
       call getenv_loc('SIDEPAR',sidename)
-      open (isidep,file=sidename,status='old',readonly)
+      open (isidep,file=sidename,status='old',action='read')
 #endif
 #ifndef OLDSCP
 C
@@ -2129,7 +2172,7 @@ C
 #elif (defined G77)
       open (iscpp,file=scpname,status='old')
 #else
-      open (iscpp,file=scpname,status='old',readonly)
+      open (iscpp,file=scpname,status='old',action='read')
 #endif
 #endif
       call getenv_loc('PATTERN',patname)
@@ -2140,7 +2183,7 @@ C
 #elif (defined G77)
       open (icbase,file=patname,status='old')
 #else
-      open (icbase,file=patname,status='old',readonly)
+      open (icbase,file=patname,status='old',action='read')
 #endif
 #ifdef MPI
 C Open output file only for CG processes
@@ -2195,6 +2238,7 @@ c      print *,"Processor",myrank," fg_rank",fg_rank
 #if defined(AIX) || defined(PGI)
       if (me.eq.king .or. .not. out1file) 
      &   open(iout,file=outname,status='unknown')
+c#define DEBUG
 #ifdef DEBUG
       if (fg_rank.gt.0) then
         write (liczba,'(i3.3)') myrank/nfgtasks
@@ -2203,6 +2247,7 @@ c      print *,"Processor",myrank," fg_rank",fg_rank
      &   status='unknown')
       endif
 #endif
+c#undef DEBUG
       if(me.eq.king) then
        open(igeom,file=intname,status='unknown',position='append')
        open(ipdb,file=pdbname,status='unknown')
@@ -2214,14 +2259,17 @@ c1out       open(iout,file=outname,status='unknown')
 #else
       if (me.eq.king .or. .not.out1file)
      &    open(iout,file=outname,status='unknown')
+c#define DEBUG
 #ifdef DEBUG
       if (fg_rank.gt.0) then
+        print "Processor",fg_rank," opening output file"
         write (liczba,'(i3.3)') myrank/nfgtasks
         write (ll,'(bz,i3.3)') fg_rank
         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
      &   status='unknown')
       endif
 #endif
+c#undef DEBUG
       if(me.eq.king) then
        open(igeom,file=intname,status='unknown',access='append')
        open(ipdb,file=pdbname,status='unknown')