Merge branch 'UCGM' of mmka.chem.univ.gda.pl:unres4 into UCGM
[unres4.git] / source / wham / io_wham.F90
index 4c8a65d..d8cb743 100644 (file)
       open (iscpp_nucl,file=scpname_nucl,status='old')
       call mygetenv('IONPAR_NUCL',ionnuclname)
       open (iionnucl,file=ionnuclname,status='old')
+      call mygetenv('IONPAR_TRAN',iontranname)
+      open (iiontran,file=iontranname,status='old',action='read')
 
 
 #ifndef OLDSCP
       call reada(controlcard,"DTRISS",dtriss,1.001D0)
       call reada(controlcard,"SSSCALE",ssscale,1.0D0)
       dyn_ss=(index(controlcard,'DYN_SS').gt.0)
-
+      call reada(controlcard,"LIPSCALE",lipscale,1.0D0)
 allocate(ww(max_eneW))
       do i=1,n_eneW
         key = wname(i)(:ilen(wname(i)))
@@ -543,6 +545,7 @@ allocate(ww(max_eneW))
       wscpho=ww(48)
       wpeppho=ww(49)
       wcatnucl=ww(50)
+      wcat_tran=ww(56)
 !      print *,"KURWA",ww(48)
 !        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   "
 !        "WVDWPP    ","WELPP     ","WVDWPSB   ","WELPSB    ","WVDWSB    ",&
@@ -600,6 +603,7 @@ allocate(ww(max_eneW))
       weights(48) =wscpho
       weights(49) =wpeppho
       weights(50) =wcatnucl
+      weights(56)=wcat_tran
 
 ! el--------
       call card_concat(controlcard,.false.)
@@ -726,8 +730,8 @@ allocate(ww(max_eneW))
           enddo
         enddo
       endif
-            if (.not. allocated(msc)) allocate(msc(ntyp1,5))
-            if (.not. allocated(restok)) allocate(restok(ntyp1,5))
+            if (.not. allocated(msc)) allocate(msc(-ntyp1:ntyp1,5))
+            if (.not. allocated(restok)) allocate(restok(-ntyp1:ntyp1,5))
        if (oldion.eq.1) then
 
             do i=1,ntyp_molec(5)
@@ -2774,7 +2778,7 @@ allocate(ww(max_eneW))
       if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
 
 
-      allocate (ichargecat(ntyp_molec(5)))
+      allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
        if (oldion.eq.0) then
             if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
@@ -2784,13 +2788,13 @@ allocate(ww(max_eneW))
              read(iion,*) ijunk
             endif
 
-            do i=1,ntyp_molec(5)
+            do i=-ntyp_molec(5),ntyp_molec(5)
              read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
              print *,msc(i,5),restok(i,5)
             enddo
             ip(5)=0.2
 !DIR$ NOUNROLL 
-      do j=1,ntyp_molec(5)
+      do j=1,ntyp_molec(5)-1
        do i=1,ntyp
 !       do j=1,ntyp_molec(5)
 !        write (*,*) "Im in ALAB", i, " ", j
@@ -2818,7 +2822,7 @@ allocate(ww(max_eneW))
       END DO
       allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
       do i=1,ntyp
-        do j=1,ntyp_molec(5)
+        do j=1,ntyp_molec(5)-1 !without zinc
           epsij=epscat(i,j)
           rrij=sigmacat(i,j)
           rrij=rrij**expon
@@ -2870,7 +2874,36 @@ allocate(ww(max_eneW))
        enddo
 
       endif
-
+! here two denotes the Zn2+ and Cu2+
+      write(iout,*) "before TRANPARM"
+      allocate(aomicattr(0:3,2))
+      allocate(athetacattran(0:6,5,2))
+      allocate(agamacattran(3,5,2))
+      allocate(acatshiftdsc(5,2))
+      allocate(bcatshiftdsc(5,2))
+      allocate(demorsecat(5,2))
+      allocate(alphamorsecat(5,2))
+      allocate(x0catleft(5,2))
+      allocate(x0catright(5,2))
+      allocate(x0cattrans(5,2))
+      allocate(ntrantyp(2))
+      do i=1,1 ! currently only Zn2+
+
+      read(iiontran,*) ntrantyp(i)
+!ntrantyp=4
+!| ao0 ao1 ao2 ao3
+!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
+!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
+!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+      read(iiontran,*) (aomicattr(j,i),j=0,3)
+      do j=1,ntrantyp(i)
+       read (iiontran,*) (agamacattran(k,j,i),k=1,3),&
+       (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+       demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+       x0cattrans(j,i)
+      enddo
+      enddo
 #ifdef CLUSTER
 !
 ! Define the SC-p interaction constants
@@ -3103,11 +3136,13 @@ allocate(ww(max_eneW))
       subroutine read_general_data(*)
 
       use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
-          scelemode,TUBEmode,tor_mode,energy_dec
+          scelemode,TUBEmode,tor_mode,energy_dec,r_cut_ang,r_cut_mart,&
+          rlamb_mart
          
       use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss,constr_homology
       use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
-          bordtubebot,tubebufthick,buftubebot,buftubetop
+          bordtubebot,tubebufthick,buftubebot,buftubetop,buflipbot, bufliptop,bordlipbot,bordliptop,     &
+        lipbufthick,lipthick
 !      implicit none
 !      include "DIMENSIONS"
 !      include "DIMENSIONS.ZSCOPT"
@@ -3126,6 +3161,7 @@ allocate(ww(max_eneW))
 !      include "COMMON.FREE"
 !      include "COMMON.CONTROL"
 !      include "COMMON.ENERGIES"
+       
       character(len=800) :: controlcard
       integer :: i,j,k,ii,n_ene_found
       integer :: ind,itype1,itype2,itypf,itypsc,itypp
@@ -3187,6 +3223,23 @@ 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)
+      call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+      call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+      if (lipthick.gt.0.0d0) then
+       bordliptop=(boxzsize+lipthick)/2.0
+       bordlipbot=bordliptop-lipthick
+      if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+      write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+      buflipbot=bordlipbot+lipbufthick
+      bufliptop=bordliptop-lipbufthick
+      if ((lipbufthick*2.0d0).gt.lipthick) &
+       write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+      endif !lipthick.gt.0
+      write(iout,*) "bordliptop=",bordliptop
+      write(iout,*) "bordlipbot=",bordlipbot
+      write(iout,*) "bufliptop=",bufliptop
+      write(iout,*) "buflipbot=",buflipbot
+
       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
       call readi(controlcard,"SCELEMODE",scelemode,0)
       call readi(controlcard,"OLDION",oldion,0)
@@ -3220,6 +3273,9 @@ allocate(ww(max_eneW))
       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
+      call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
+      call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
+      call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
       check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
       call readi(controlcard,'SYM',symetr,1)
@@ -3972,8 +4028,8 @@ allocate(ww(max_eneW))
       do i=nnt,nct
         mnum=molnum(i)
         iti=itype(i,mnum)
-        if (iti.eq.ntyp1) then
-          if (itype(i-1,molnum(i-1)).eq.ntyp1) then
+        if (iti.eq.ntyp1_molec(mnum)) then
+          if (itype(i-1,molnum(i-1)).eq.ntyp1_molec(mnum)) then
           ichain=ichain+1
           ires=0
           write (ipdb,'(a)') 'TER'
@@ -4002,7 +4058,7 @@ allocate(ww(max_eneW))
       write (ipdb,'(a)') 'TER'
       do i=nnt,nct-1
         mnum=molnum(i)
-        if (itype(i,mnum).eq.ntyp1) cycle
+        if (itype(i,mnum).eq.ntyp1_molec(mnum)) cycle
         if (itype(i,mnum).eq.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then
           write (ipdb,30) ica(i),ica(i+1)
         else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum).ne.ntyp1_molec(mnum)) then