Calculation of folded fractions (at given RMSD cutoff) in WHAM analysis (7th column...
authorPawel Krupa <vetinari@piasek4.chem.univ.gda.pl>
Thu, 10 Mar 2016 22:32:05 +0000 (23:32 +0100)
committerPawel Krupa <vetinari@piasek4.chem.univ.gda.pl>
Thu, 10 Mar 2016 22:32:05 +0000 (23:32 +0100)
source/wham/src-M/COMMON.ENERGIES
source/wham/src-M/COMMON.PROTFILES
source/wham/src-M/DIMENSIONS.FREE
source/wham/src-M/enecalc1.F
source/wham/src-M/readrtns.F
source/wham/src-M/wham_calc1.F

index 2d40a95..1c82395 100644 (file)
@@ -1,4 +1,4 @@
       double precision potE(MaxStr_Proc,Max_Parm),entfac(MaxStr_Proc),
-     & q(MaxQ+2,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm)
+     & q(MaxQ+3,MaxStr_Proc),enetb(max_ene,MaxStr_Proc,Max_Parm)
       integer einicheck
       common /energies/ potE,entfac,q,enetb,einicheck
index 3287326..f5a1b80 100644 (file)
@@ -5,6 +5,7 @@
      &  nfile_cx(MaxR,MaxT_h,Max_Parm),
      &  rec_start(MaxR,MaxT_h,Max_Parm),
      &  rec_end(MaxR,MaxT_h,Max_Parm),lenrec,lenrec1,lenrec2
-      common /protfil/ protfiles,bprotfiles,
+      real*8 frac_cutoff
+      common /protfil/ protfiles,bprotfiles, frac_cutoff,
      &  nfile_bin,nfile_asc,nfile_cx,rec_start,rec_end,lenrec,lenrec1,
      &  lenrec2
index 691d9b2..304e83d 100644 (file)
@@ -3,7 +3,7 @@
       integer MaxR,MaxT_h
       integer MaxSlice
       parameter (Max_Parm=1)
-      parameter (MaxQ=1,MaxQ1=MaxQ+2)
+      parameter (MaxQ=1,MaxQ1=MaxQ+3)
       parameter(MaxR=1,MaxT_h=32)
       parameter(MaxSlice=40)
       integer MaxN
index c4f1f45..697c31d 100644 (file)
            q(nQ+1,iii+1)=rmsnat(iii+1)
          endif
          q(nQ+2,iii+1)=gyrate(iii+1)
+         if (q(nQ+1,iii+1).le.frac_cutoff) then
+           q(nQ+3,iii+1)=1.0d0
+         else
+           q(nQ+3,iii+1)=0.0d0
+         endif
 c        fT=T0*beta_h(ib,ipar)*1.987D-3
 c        ft=2.0d0/(1.0d0+1.0d0/(T0*beta_h(ib,ipar)*1.987D-3))
         if (rescale_mode.eq.1) then
index 5dd092e..6762578 100644 (file)
@@ -51,6 +51,8 @@
       call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
       write (iout,*) "MaxSlice",MaxSlice
       call readi(controlcard,"NSLICE",nslice,1)
+      call reada(controlcard,"FRAC_CUTOFF",frac_cutoff,10.0d0)
+      write (iout,*) "RMSD cutoff for fraction calculations",frac_cutoff
       call flush(iout)
       if (nslice.gt.MaxSlice) then
         write (iout,*) "Error: parameter out of range: NSLICE",nslice,
index 6768dcf..ff78e8d 100644 (file)
@@ -531,7 +531,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           sumE_p(i,iparm)=0.0d0
           sumEbis_p(i,iparm)=0.0d0
           sumEsq_p(i,iparm)=0.0d0
-          do j=1,nQ+2
+          do j=1,nQ+3
             sumQ_p(j,i,iparm)=0.0d0
             sumQsq_p(j,i,iparm)=0.0d0
             sumEQ_p(j,i,iparm)=0.0d0
@@ -546,7 +546,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           sumE(i,iparm)=0.0d0
           sumEbis(i,iparm)=0.0d0
           sumEsq(i,iparm)=0.0d0
-          do j=1,nQ+2
+          do j=1,nQ+3
             sumQ(j,i,iparm)=0.0d0
             sumQsq(j,i,iparm)=0.0d0
             sumEQ(j,i,iparm)=0.0d0
@@ -826,7 +826,7 @@ c            write (iout,*) "ftbis",ftbis
             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
-            do j=1,nQ+2
+            do j=1,nQ+3
               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
@@ -837,7 +837,7 @@ c            write (iout,*) "ftbis",ftbis
             sumE(k,iparm)=sumE(k,iparm)+etot*weight
             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
-            do j=1,nQ+2
+            do j=1,nQ+3
               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
@@ -1087,7 +1087,7 @@ c            write (iout,*) "ftbis",ftbis
      &    sumW(i,iparm)
         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
-        do j=1,nQ+2
+        do j=1,nQ+3
           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
      &     -sumQ(j,i,iparm)**2
@@ -1098,15 +1098,15 @@ c            write (iout,*) "ftbis",ftbis
      &   (startGridT+i*delta_T))+potEmin
         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
      &   sumW(i,iparm),sumE(i,iparm)
-        write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+        write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+3)
         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
-     &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+     &   (sumQsq(j,i,iparm),j=1,nQ+3),(sumEQ(j,i,iparm),j=1,nQ+3)
         write (iout,*) 
         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
      &   sumW(i,iparm),sumE(i,iparm)
-        write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+        write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+3)
         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
-     &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+     &   (sumQsq(j,i,iparm),j=1,nQ+3),(sumEQ(j,i,iparm),j=1,nQ+3)
         write (34,*) 
         call flush(34)
       enddo