Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / wham / io_wham.F90
index 97b74d1..fee1e46 100644 (file)
@@ -98,6 +98,8 @@
 
       call mygetenv('SCPPAR_NUCL',scpname_nucl)
       open (iscpp_nucl,file=scpname_nucl,status='old')
+      call mygetenv('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old')
 
 
 #ifndef OLDSCP
 !
       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)
           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
       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
       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)))
@@ -485,6 +516,7 @@ allocate(ww(max_eneW))
       wpepbase=ww(47)
       wscpho=ww(48)
       wpeppho=ww(49)
+      wcatnucl=ww(50)
 !      print *,"KURWA",ww(48)
 !        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   "
 !        "WVDWPP    ","WELPP     ","WVDWPSB   ","WELPSB    ","WVDWSB    ",&
@@ -541,6 +573,8 @@ allocate(ww(max_eneW))
       weights(47)= wpepbase
       weights(48) =wscpho
       weights(49) =wpeppho
+      weights(50) =wcatnucl
+
 ! el--------
       call card_concat(controlcard,.false.)
 
@@ -584,7 +618,7 @@ allocate(ww(max_eneW))
       write (iout,*) "Parameter set:",iparm
       write (iout,*) "Energy-term weights:"
       do i=1,n_eneW
-        write (iout,'(a16,f10.5)') wname(i),ww(i)
+        write (iout,'(i3,a16,f10.5)') i,wname(i),ww(i)
       enddo
       write (iout,*) "Sidechain potential file        : ",&
         sidename_t(:ilen(sidename_t))
@@ -683,6 +717,23 @@ allocate(ww(max_eneW))
             enddo
             print *, catprm
       endif
+      allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
+      do i=1,ntyp_molec(5)
+         do j=1,ntyp_molec(2)
+         write(iout,*) i,j
+            read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
+         enddo
+      enddo
+      write(*,'(3(5x,a6)11(7x,a6))') "w1    ","w2    ","epslj ","pis1  ", &
+      "sigma0","epsi0 ","chi1   ","chip1 ","sig   ","b1    ","b2    ", &
+      "b3    ","b4    ","chis1  "
+      do i=1,ntyp_molec(5)
+         do j=1,ntyp_molec(2)
+            write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
+                                      restyp(i,5),"-",restyp(j,2)
+         enddo
+      enddo
+
       read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
       do i=1,ntyp_molec(2)
         nbondterm_nucl(i)=1
@@ -2677,6 +2728,7 @@ allocate(ww(max_eneW))
 ! Ions by Aga
 
        allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+       allocate(alphapolcat2(ntyp,ntyp))
        allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
        allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
        allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
@@ -2733,7 +2785,7 @@ allocate(ww(max_eneW))
 !       rborncat(i,j),rborncat(j,i),&
        rborn1cat(i,j),rborn2cat(i,j),&
        (wqdipcat(k,i,j),k=1,2), &
-       alphapolcat(i,j),alphapolcat(j,i), &
+       alphapolcat(i,j),alphapolcat2(j,i), &
        (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
 !       print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) 
        END DO
@@ -2913,7 +2965,7 @@ allocate(ww(max_eneW))
 !
 ! Define the constants of the disulfide bridge
 !
-      ebr=-5.50D0
+!      ebr=-5.50D0
 !
 ! Old arbitrary potential - commented out.
 !
@@ -2924,14 +2976,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
@@ -3013,7 +3079,7 @@ allocate(ww(max_eneW))
       use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
           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
@@ -3098,7 +3164,7 @@ allocate(ww(max_eneW))
       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
@@ -3842,8 +3908,9 @@ allocate(ww(max_eneW))
 !--------------------------------------------------------------------------------
       subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
 
-      use geometry_data, only:nres,c
+      use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize
       use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
+      use energy, only:boxshift
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 !      include 'DIMENSIONS.ZSCOPT'
@@ -3858,7 +3925,7 @@ allocate(ww(max_eneW))
                       'D','E','F','G','H','I','J','K','L','M','N','O',&
                       'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
       integer,dimension(nres) :: ica !(maxres)
-      real(kind=8) :: temp,efree,etot,entropy,rmsdev
+      real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj
       integer :: ii,i,j,iti,ires,iatom,ichain,mnum
       write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
         ii,temp,rmsdev
@@ -3882,9 +3949,17 @@ allocate(ww(max_eneW))
         ires=ires+1
         iatom=iatom+1
         ica(i)=iatom
+        if (mnum.ne.5) then
         write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
            ires,(c(j,i),j=1,3)
-        if (iti.ne.10) then
+        else
+        xj=boxshift(c(1,i)-c(1,2),boxxsize)
+        yj=boxshift(c(2,i)-c(2,2),boxysize)
+        zj=boxshift(c(3,i)-c(3,2),boxzsize)
+        write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
+           ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj
+        endif
+        if ((iti.ne.10).and.(mnum.ne.5)) then
           iatom=iatom+1
           write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
             ires,(c(j,nres+i),j=1,3)