correcation in ion param and dyn_ss

[unres4.git] / source / wham / io_wham.F90

diff --git a/source/wham/io_wham.F90 b/source/wham/io_wham.F90
index b96b816..e3f0151 100644 (file)
--- a/source/wham/io_wham.F90
+++ b/source/wham/io_wham.F90
@@ -163,9 +163,9 @@
 !
       use energy_data
       use geometry_data, only:nres,deg2rad,c,dc,nres_molec
-      use control_data, only:iscode
+      use control_data, only:iscode,dyn_ss
       use io_base, only:rescode
-      use control, only:setup_var,init_int_table
+      use control, only:setup_var,init_int_table,hpb_partition
       use geometry, only:alloc_geo_arrays
       use energy, only:alloc_ener_arrays      
 !      implicit real*8 (a-h,o-z)
@@ -344,6 +344,20 @@
           dhpb(i),ebr,forcon(i)
         enddo
       endif
+      if (ns.gt.0.and.dyn_ss) then
+          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
+      endif
       write (iout,'(a)')
       return
       end subroutine molread
@@ -410,14 +424,14 @@
       real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
                 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,    &
                 epspeptube,epssctube,sigmapeptube,      &
-                sigmasctube
+                sigmasctube,ssscale
 
       real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
                 res1
       integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
       integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
       character*3 string
-
+      
 !
 ! Body
 !
@@ -429,10 +443,26 @@
       vbl=3.8D0
       vblinv=1.0D0/vbl
       vblinv2=vblinv*vblinv
+      itime_mat=0
 #ifndef CLUSTER
       call card_concat(controlcard,.true.)
       wname(4)="WCORRH"
 !el
+      call reada(controlcard,"D0CM",d0cm,3.78d0)
+      call reada(controlcard,"AKCM",akcm,15.1d0)
+      call reada(controlcard,"AKTH",akth,11.0d0)
+      call reada(controlcard,"AKCT",akct,12.0d0)
+      call reada(controlcard,"V1SS",v1ss,-1.08d0)
+      call reada(controlcard,"V2SS",v2ss,7.61d0)
+      call reada(controlcard,"V3SS",v3ss,13.7d0)
+      call reada(controlcard,"EBR",ebr,-5.50D0)
+      call reada(controlcard,"ATRISS",atriss,0.301D0)
+      call reada(controlcard,"BTRISS",btriss,0.021D0)
+      call reada(controlcard,"CTRISS",ctriss,1.001D0)
+      call reada(controlcard,"DTRISS",dtriss,1.001D0)
+      call reada(controlcard,"SSSCALE",ssscale,1.0D0)
+      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+
 allocate(ww(max_eneW))
       do i=1,n_eneW
         key = wname(i)(:ilen(wname(i)))
@@ -2912,7 +2942,7 @@ allocate(ww(max_eneW))
 !
 ! Define the constants of the disulfide bridge
 !
-      ebr=-5.50D0
+!      ebr=-5.50D0
 !
 ! Old arbitrary potential - commented out.
 !
@@ -2923,14 +2953,28 @@ allocate(ww(max_eneW))
 ! energy surface of diethyl disulfide.
 ! A. Liwo and U. Kozlowska, 11/24/03
 !
-      D0CM = 3.78d0
-      AKCM = 15.1d0
-      AKTH = 11.0d0
-      AKCT = 12.0d0
-      V1SS =-1.08d0
-      V2SS = 7.61d0
-      V3SS = 13.7d0
+!      D0CM = 3.78d0
+!      AKCM = 15.1d0
+!      AKTH = 11.0d0
+!      AKCT = 12.0d0
+!      V1SS =-1.08d0
+!      V2SS = 7.61d0
+!      V3SS = 13.7d0
+#ifndef CLUSTER
 
+      if (dyn_ss) then
+       ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
+        Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
+        akcm=akcm/wsc*ssscale
+        akth=akth/wsc*ssscale
+        akct=akct/wsc*ssscale
+        v1ss=v1ss/wsc*ssscale
+        v2ss=v2ss/wsc*ssscale
+        v3ss=v3ss/wsc*ssscale
+      else
+        ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+      endif
+#endif
       if (lprint) then
       write (iout,'(/a)') "Disulfide bridge parameters:"
       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
@@ -3010,9 +3054,9 @@ allocate(ww(max_eneW))
       subroutine read_general_data(*)
 
       use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
-          scelemode,TUBEmode,tor_mode
+          scelemode,TUBEmode,tor_mode,energy_dec
          
-      use energy_data, only:distchainmax,tubeR0,tubecenter
+      use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss
       use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
           bordtubebot,tubebufthick,buftubebot,buftubetop
 !      implicit none
@@ -3094,9 +3138,10 @@ allocate(ww(max_eneW))
       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+      energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
       call readi(controlcard,"SCELEMODE",scelemode,0)
       call readi(controlcard,"OLDION",oldion,0)
-
+      dyn_ss=(index(controlcard,'DYN_SS').gt.0)
       print *,"SCELE",scelemode
       call readi(controlcard,'TORMODE',tor_mode,0)
 !C      if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
@@ -3116,7 +3161,7 @@ allocate(ww(max_eneW))
        buftubetop=bordtubetop-tubebufthick
       endif
       ions=index(controlcard,"IONS").gt.0
-      call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+      call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
       call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
       write(iout,*) "R_CUT_ELE=",r_cut_ele
       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0