Implementation of ions into NARES

diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt
index 3b6f103..8852d5f 100644 (file)
--- a/source/unres/CMakeLists.txt
+++ b/source/unres/CMakeLists.txt
@@ -71,7 +71,7 @@ set(UNRES_MD_SRC3
 if (Fortran_COMPILER_NAME STREQUAL "ifort")
   set (CMAKE_Fortran_FLAGS_RELEASE " ")
   set (CMAKE_Fortran_FLAGS_DEBUG   "-O0 -g -traceback")
-#  set(FFLAGS0 "-fpp -c -O3 -ip " ) 
+#  set(FFLAGS0 "-fpp -c -CB -g -ip " ) 
   set(FFLAGS0 "-O3 -ip -fpp  -heap-arrays -mcmodel=medium" ) 
 #  set(FFLAGS0 "-O0 -CB -CA -g" )
   set(FFLAGS1 "-fpp -c -O " ) 

diff --git a/source/unres/control.F90 b/source/unres/control.F90
index b46dd40..1952c9a 100644 (file)
--- a/source/unres/control.F90
+++ b/source/unres/control.F90
@@ -247,6 +247,7 @@
       itube=61
 !     IONS
       iion=401
+      iionnucl=402
 #if defined(WHAM_RUN) || defined(CLUSTER)
 !
 ! setting the mpi variables for WHAM

diff --git a/source/unres/data/energy_data.F90 b/source/unres/data/energy_data.F90
index 8fd85cd..fc3ba37 100644 (file)
--- a/source/unres/data/energy_data.F90
+++ b/source/unres/data/energy_data.F90
@@ -73,7 +73,7 @@
        wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, &
        wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,&
        wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb,wcatprot,&
-       wcatcat,wscbase,wpepbase,wscpho,wpeppho,wdihc
+       wcatcat,wscbase,wpepbase,wscpho,wpeppho,wdihc,wcatnucl
 #ifdef CLUSTER
       real(kind=8) :: scalscp
 #endif
@@ -364,6 +364,8 @@
        real(kind=8) :: distafminit,forceAFMconst,velAFMconst
       integer :: afmend,afmbeg
       real(kind=8),dimension(:,:), allocatable :: catprm
+      real(kind=8),dimension(:,:,:), allocatable :: catnuclprm
+
 
          real(kind=8),dimension(:,:), allocatable ::  eps_scbase, &
         sigma_scbase,                         &

diff --git a/source/unres/data/io_units.F90 b/source/unres/data/io_units.F90
index 909e6fd..764889b 100644 (file)
--- a/source/unres/data/io_units.F90
+++ b/source/unres/data/io_units.F90
@@ -17,7 +17,7 @@
        ithep_pdb,irotam_pdb,iliptranpar,itube,   &
        ibond_nucl,ithep_nucl,irotam_nucl,itorp_nucl,itordp_nucl,     &
        isidep_nucl,iscpp_nucl,isidep_scbase,isidep_pepbase,isidep_scpho,&
-       isidep_peppho,iion      
+       isidep_peppho,iion,iionnucl      
 #ifdef WHAM_RUN
 ! el wham iounits
       integer :: isidep1,ihist,iweight,izsc,idistr
@@ -50,7 +50,7 @@
        bondname_nucl,thetname_nucl,rotname_nucl,torname_nucl, &
        tordname_nucl,sidename_nucl,scpname_nucl,              &
        sidename_scbase,pepname_pepbase,pepname_scpho,pepname_peppho, &
-       ionname
+       ionname,ionnuclname
 !-----------------------------------------------------------------------
 ! INP    - main input file
 ! IOUT   - list file

diff --git a/source/unres/data/names.F90 b/source/unres/data/names.F90
index bb5a975..d1b8d11 100644 (file)
--- a/source/unres/data/names.F90
+++ b/source/unres/data/names.F90
@@ -66,7 +66,7 @@
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
 ! Number of energy components
-      integer,parameter :: n_ene=49
+      integer,parameter :: n_ene=50
       integer :: n_ene2=2*n_ene
 !-----------------------------------------------------------------------------
 ! common.names
@@ -95,7 +95,8 @@
         "EESSB     ","ESTR      ","EBE       ","ESBLOC    ","ETORS     ",&
         "ETORSD    ","ECORR     ","ECORR3    ","NULL      ","NULL      ",&
         "ECATPROT  ","ECATCAT   ","NULL      ","NULL      ","NULL      ",&
-        "ESCBASE   ","EPEPBASE  ","ESCPHO    ","EPEPPHO   "/)
+        "ESCBASE   ","EPEPBASE  ","ESCPHO    ","EPEPPHO   ",&
+        "ECATION_NUCL"/)
 
       character(len=10),dimension(n_ene) :: wname = &
       (/"WSC       ","WSCP      ","WELEC     ","WCORR     ","WCORR5    ","WCORR6    ","WEL_LOC   ",&
@@ -106,13 +107,13 @@
         "WELSB     ","WBOND_NUCL","WANG_NUCL ","WSBLOC    ","WTOR_NUCL ",&
         "WTORD_NUCL","WCORR_NUCL","WCORR3_NUC","WNULL     ","WNULL     ",&
         "WCATPROT  ","WCATCAT   ","WNULL     ","WNULL     ","WNULL     ",&
-        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   "/)
+        "WSCBASE   ","WPEPBASE  ","WSCPHO    ","WPEPPHO   ","WCATNUCL  "/)
 
       integer :: nprint_ene = 21
       integer,dimension(n_ene) :: print_order = &
          (/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25,&
          26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,&
-         48,49/)
+         48,49,50/)
 
       character(len=1), dimension(2) :: sugartyp = (/'D',' '/)
 !#endif

diff --git a/source/unres/energy.F90 b/source/unres/energy.F90
index 350a50d..376d318 100644 (file)
--- a/source/unres/energy.F90
+++ b/source/unres/energy.F90
@@ -137,7 +137,7 @@
          gvdwc_peppho
 !------------------------------IONS GRADIENT
         real(kind=8),dimension(:,:),allocatable  ::  gradcatcat, &
-          gradpepcat,gradpepcatx
+          gradpepcat,gradpepcatx,gradnuclcat,gradnuclcatx
 !      real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2)
 
 
@@ -247,7 +247,8 @@
                       ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
                       ecorr3_nucl
 ! energies for ions 
-      real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
+      real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
+                      ecation_nucl
 ! energies for protein nucleic acid interaction
       real(kind=8) :: escbase,epepbase,escpho,epeppho
 
@@ -333,6 +334,7 @@
           weights_(47)=wpepbase
           weights_(48)=wscpho
           weights_(49)=wpeppho
+          weights_(50)=wcatnucl          
 !          wcatcat= weights(41)
 !          wcatprot=weights(42)
 
@@ -381,6 +383,7 @@
           wpepbase=weights(47)
           wscpho=weights(48)
           wpeppho=weights(49)
+          wcatnucl=weights(50)
 !      welpsb=weights(28)*fact(1)
 !
 !      wcorr_nucl= weights(37)*fact(1)
@@ -398,12 +401,12 @@
         call MPI_Bcast(itime_mat,1,MPI_INT,king,FG_COMM,IERROR)
         endif
        if (mod(itime_mat,imatupdate).eq.0) call make_SCp_inter_list
-       write (iout,*) "after make_SCp_inter_list"
+!       write (iout,*) "after make_SCp_inter_list"
        if (mod(itime_mat,imatupdate).eq.0) call make_SCSC_inter_list
-       write (iout,*) "after make_SCSC_inter_list"
+!       write (iout,*) "after make_SCSC_inter_list"
 
        if (mod(itime_mat,imatupdate).eq.0) call make_pp_inter_list
-       write (iout,*) "after make_pp_inter_list"
+!       write (iout,*) "after make_pp_inter_list"
 
 !      print *,'Processor',myrank,' calling etotal ipot=',ipot
 !      print *,'Processor',myrank,' nnt=',nnt,' nct=',nct
@@ -841,6 +844,7 @@
       call epsb(evdwpsb,eelpsb)
       call esb(esbloc)
       call multibody_hb_nucl(ecorr_nucl,ecorr3_nucl,n_corr,n_corr1)
+            call ecat_nucl(ecation_nucl)
       else
        etors_nucl=0.0d0
        estr_nucl=0.0d0
@@ -855,6 +859,7 @@
        evdwpp=0.0d0
        eespp=0.0d0
        etors_d_nucl=0.0d0
+       ecation_nucl=0.0d0
       endif
 !      write(iout,*) ecorr_nucl,"ecorr_nucl",nres_molec(2)
 !      print *,"before ecatcat",wcatcat
@@ -958,6 +963,7 @@
       energia(48)=escpho
       energia(49)=epeppho
 !      energia(50)=ecations_prot_amber
+      energia(50)=ecation_nucl
       call sum_energy(energia,.true.)
       if (dyn_ss) call dyn_set_nss
 !      print *," Processor",myrank," left SUM_ENERGY"
@@ -1000,7 +1006,8 @@
       real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
                       ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
                       ecorr3_nucl
-      real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
+      real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
+                      ecation_nucl
       real(kind=8) :: escbase,epepbase,escpho,epeppho
       integer :: i
 #ifdef MPI
@@ -1084,6 +1091,7 @@
       epepbase=energia(47)
       escpho=energia(48)
       epeppho=energia(49)
+      ecation_nucl=energia(50)
 !      ecations_prot_amber=energia(50)
 
 !      energia(41)=ecation_prot
@@ -1103,7 +1111,7 @@
        +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
        +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
        +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
-       +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+       +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) &
        +wang*ebe+wtor*etors+wscloc*escloc &
@@ -1117,7 +1125,7 @@
        +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
        +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
        +wcatprot*ecation_prot+wcatcat*ecationcation+wscbase*escbase&
-       +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho
+       +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl
 #endif
       energia(0)=etot
 ! detecting NaNQ
@@ -1250,7 +1258,8 @@
       real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,&
                       ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,&
                       ecorr3_nucl
-      real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber
+      real(kind=8) :: ecation_prot,ecationcation,ecations_prot_amber,&
+                      ecation_nucl
       real(kind=8) :: escbase,epepbase,escpho,epeppho
 
       etot=energia(0)
@@ -1304,6 +1313,7 @@
       epepbase=energia(47)
       escpho=energia(48)
       epeppho=energia(49)
+      ecation_nucl=energia(50)
 !      ecations_prot_amber=energia(50)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,&
@@ -1320,7 +1330,7 @@
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
         ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat, &
         escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
-        etot
+        ecation_nucl,wcatnucl,etot
    10 format (/'Virtual-chain energies:'// &
        'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
        'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
@@ -1367,6 +1377,7 @@
        'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ &
        'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/&
        'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/&
+       'ECATBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(cation nucl-base)'/&
        'ETOT=  ',1pE16.6,' (total)')
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,&
@@ -1383,7 +1394,7 @@
         etors_d_nucl,wtor_d_nucl,ecorr_nucl,wcorr_nucl,&
         ecorr3_nucl,wcorr3_nucl,ecation_prot,wcatprot,ecationcation,wcatcat,  &
         escbase,wscbase,epepbase,wpepbase,escpho,wscpho,epeppho,wpeppho,&
-        etot
+        ecation_nucl,wcatnucl,etot
    10 format (/'Virtual-chain energies:'// &
        'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ &
        'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ &
@@ -1429,6 +1440,7 @@
        'EPEPBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-base)'/ &
        'ESCPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(sc-prot nucl-phosphate)'/&
        'EPEPPHO=',1pE16.6,' WEIGHT=',1pD16.6,'(pep-prot nucl-phosphate)'/&
+       'ECATBASE=',1pE16.6,' WEIGHT=',1pD16.6,'(cation nucl-base)'/&
        'ETOT=  ',1pE16.6,' (total)')
 #endif
       return
@@ -1902,12 +1914,12 @@
       dCAVdOM1=0.0d0 
       dGCLdOM1=0.0d0 
       dPOLdOM1=0.0d0
-
+!             write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j
 
       do icont=g_listscsc_start,g_listscsc_end
       i=newcontlisti(icont)
       j=newcontlistj(icont)
-
+!      write (iout,*) "RWA", g_listscsc_start,g_listscsc_end,i,j
 !      do i=iatsc_s,iatsc_e
 !C        print *,"I am in EVDW",i
         itypi=iabs(itype(i,1))
@@ -1995,6 +2007,7 @@
            zj=c(3,nres+j)
               call to_box(xj,yj,zj)
               call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+              write (iout,*) "KWA2", itypi,itypj
               aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
                +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
               bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
@@ -11143,7 +11156,7 @@
                      wscbase*gvdwc_scbase(j,i)+ &
                      wpepbase*gvdwc_pepbase(j,i)+&
                      wscpho*gvdwc_scpho(j,i)+   &
-                     wpeppho*gvdwc_peppho(j,i)
+                     wpeppho*gvdwc_peppho(j,i)+wcatnucl*gradnuclcat(j,i)
 
        
 
@@ -11181,7 +11194,7 @@
                      wscbase*gvdwc_scbase(j,i)+ &
                      wpepbase*gvdwc_pepbase(j,i)+&
                      wscpho*gvdwc_scpho(j,i)+&
-                     wpeppho*gvdwc_peppho(j,i)
+                     wpeppho*gvdwc_peppho(j,i)+wcatnucl*gradnuclcat(j,i)
 
 
         enddo
@@ -11422,7 +11435,7 @@
                      +wbond_nucl*gradb_nucl(j,i) &
                      +0.5d0*(wvdwpp_nucl*gvdwpp_nucl(j,i)+welpp*gelpp(j,i)&
                      +wvdwpsb*gvdwpsb1(j,i))&
-                     +wsbloc*gsbloc(j,i)
+                     +wsbloc*gsbloc(j,i)+wcatnucl*gradnuclcat(j,i)
 
 
 
@@ -11449,7 +11462,7 @@
                        +wcatprot* gradpepcatx(j,i)&
                        +wscbase*gvdwx_scbase(j,i) &
                        +wpepbase*gvdwx_pepbase(j,i)&
-                       +wscpho*gvdwx_scpho(j,i)
+                       +wscpho*gvdwx_scpho(j,i)+wcatnucl*gradnuclcatx(j,i)
 !              if (i.eq.3) print *,"tu?", wscpho,gvdwx_scpho(j,i)
 
         enddo
@@ -16924,6 +16937,8 @@ chip1=chip(itypi)
             gvdwx_scpho(j,i)=0.0d0
             gvdwc_scpho(j,i)=0.0d0
             gvdwc_peppho(j,i)=0.0d0
+            gradnuclcatx(j,i)=0.0d0
+            gradnuclcat(j,i)=0.0d0
           enddo
            enddo
           do i=0,nres
@@ -20343,6 +20358,8 @@ chip1=chip(itypi)
       allocate(gradpepcat(3,-1:nres))
       allocate(gradpepcatx(3,-1:nres))
       allocate(gradcatcat(3,-1:nres))
+      allocate(gradnuclcat(3,-1:nres))
+      allocate(gradnuclcatx(3,-1:nres))
 !(3,maxres)
       allocate(grad_shield_side(3,maxcontsshi,-1:nres))
       allocate(grad_shield_loc(3,maxcontsshi,-1:nres))
@@ -21155,11 +21172,11 @@ chip1=chip(itypi)
          yj=c(2,nres+j)
          zj=c(3,nres+j)
      call to_box(xj,yj,zj)
-     call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!     call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
@@ -22182,11 +22199,11 @@ chip1=chip(itypi)
          yj=c(2,j)
          zj=c(3,j)
       call to_box(xj,yj,zj)
-      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
@@ -22296,11 +22313,11 @@ chip1=chip(itypi)
          zj=c(3,j)
  
       call to_box(xj,yj,zj)
-      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
@@ -22533,7 +22550,7 @@ chip1=chip(itypi)
           Qj=Qj*2
           Qij=Qij*2
          endif
-         write(iout,*) "KURWA0",d1
+!         write(iout,*) "KURWA0",d1
 
          CALL edq_cat(ecl, elj, epol)
         eheadtail = ECL + elj + epol
@@ -23396,6 +23413,214 @@ chip1=chip(itypi)
        end subroutine ecat_prot
 
 !----------------------------------------------------------------------------
+!---------------------------------------------------------------------------
+       subroutine ecat_nucl(ecation_nucl)
+       integer i,j,k,subchap,itmp,inum,itypi,itypj
+       real(kind=8) :: xi,yi,zi,xj,yj,zj
+       real(kind=8) xj_temp,yj_temp,zj_temp,xj_safe,yj_safe,zj_safe, &
+       dist_init,dist_temp,ecation_nucl,Evan1,Evan2,Ecav,Egb,wdip1,wdip2, &
+       wvan1,wvan2,wgbsig,wgbeps,wgbchi,wgbchip,wcav1,wcav2,wcav3,wcav4, &
+       wcavsig,wcavchi,v1m,v1dpdx,wh2o,wc,Edip,rcs2,invrcs6,invrcs8,invrcs12, &
+       invrcs14,rcb,rcb2,invrcb,invrcb2,invrcb4,invrcb6,cosinus,cos2,dcosdcatconst, &
+       dcosdcalpconst,dcosdcmconst,rcav,rcav11,rcav12,constcav1,constcav2, &
+       constgb1,constgb2,constdvan1,constdvan2,sgb,sgb6,sgb7,sgb12,sgb13, &
+       cavnum,cavdenom,invcavdenom2,dcavnumdcos,dcavnumdr,dcavdenomdcos, &
+       dcavdenomdr,sslipi,ssgradlipi,sslipj,ssgradlipj,aa,bb
+       real(kind=8),dimension(3) ::gg,r,dEtotalCm,dEtotalCalp,dEvan1Cm,&
+       dEvan2Cm,cm1,cm,vcat,vsug,v1,v2,dx,vcm,dEdipCm,dEdipCalp, &
+       dEvan1Calp,dEvan2Cat,dEvan2Calp,dEtotalCat,dEdipCat,dEvan1Cat,dcosdcat, &
+       dcosdcalp,dcosdcm,dEgbdCat,dEgbdCalp,dEgbdCm,dEcavdCat,dEcavdCalp, &
+       dEcavdCm
+       real(kind=8),dimension(14) :: vcatnuclprm
+       ecation_nucl=0.0d0
+       if (nres_molec(5).eq.0) return
+       itmp=0
+       do i=1,4
+          itmp=itmp+nres_molec(i)
+       enddo
+       do i=iatsc_s_nucl,iatsc_e_nucl
+          if ((itype(i,2).eq.ntyp1_molec(2))) cycle ! leave dummy atoms
+          xi=(c(1,i+nres))
+          yi=(c(2,i+nres))
+          zi=(c(3,i+nres))
+      call to_box(xi,yi,zi)
+      call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
+          do k=1,3
+             cm1(k)=dc(k,i+nres)
+          enddo
+          do j=itmp+1,itmp+nres_molec(5)
+             xj=c(1,j)
+             yj=c(2,j)
+             zj=c(3,j)
+      call to_box(xj,yj,zj)
+!      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+      xj=boxshift(xj-xi,boxxsize)
+      yj=boxshift(yj-yi,boxysize)
+      zj=boxshift(zj-zi,boxzsize)
+
+             dist_init=xj**2+yj**2+zj**2
+
+             itypi=itype(i,2)
+             itypj=itype(j,5)
+             do k=1,13
+                vcatnuclprm(k)=catnuclprm(k,itypi,itypj)
+             enddo
+             do k=1,3
+                vcm(k)=c(k,i+nres)
+                vsug(k)=c(k,i)
+                vcat(k)=c(k,j)
+             enddo
+             do k=1,3
+                dx(k) = vcat(k)-vcm(k)
+             enddo
+             do k=1,3
+                v1(k)=dc(k,i+nres)
+                v2(k)=(vcat(k)-vsug(k))
+             enddo
+             v1m = sqrt(v1(1)**2+v1(2)**2+v1(3)**2)
+             v1dpdx = v1(1)*dx(1)+v1(2)*dx(2)+v1(3)*dx(3)
+!  The weights of the energy function calculated from
+!The quantum mechanical Gaussian simulations of potassium and sodium with deoxynucleosides
+             wh2o=78
+             wdip1 = vcatnuclprm(1)
+             wdip1 = wdip1/wh2o                     !w1
+             wdip2 = vcatnuclprm(2)
+             wdip2 = wdip2/wh2o                     !w2
+             wvan1 = vcatnuclprm(3)
+             wvan2 = vcatnuclprm(4)                 !pis1
+             wgbsig = vcatnuclprm(5)                !sigma0
+             wgbeps = vcatnuclprm(6)                !epsi0
+             wgbchi = vcatnuclprm(7)                !chi1
+             wgbchip = vcatnuclprm(8)               !chip1
+             wcavsig = vcatnuclprm(9)               !sig
+             wcav1 = vcatnuclprm(10)                !b1
+             wcav2 = vcatnuclprm(11)                !b2
+             wcav3 = vcatnuclprm(12)                !b3
+             wcav4 = vcatnuclprm(13)                !b4
+             wcavchi = vcatnuclprm(14)              !chis1
+             rcs2 = v2(1)**2+v2(2)**2+v2(3)**2
+             invrcs6 = 1/rcs2**3
+             invrcs8 = invrcs6/rcs2
+             invrcs12 = invrcs6**2
+             invrcs14 = invrcs12/rcs2
+             rcb2 = dx(1)**2+dx(2)**2+dx(3)**2
+             rcb = sqrt(rcb2)
+             invrcb = 1/rcb
+             invrcb2 = invrcb**2
+             invrcb4 = invrcb2**2
+             invrcb6 = invrcb4*invrcb2
+             cosinus = v1dpdx/(v1m*rcb)
+             cos2 = cosinus**2
+             dcosdcatconst = invrcb2/v1m
+             dcosdcalpconst = invrcb/v1m**2
+             dcosdcmconst = invrcb2/v1m**2
+             do k=1,3
+                dcosdcat(k) = (v1(k)*rcb-dx(k)*v1m*cosinus)*dcosdcatconst
+                dcosdcalp(k) = (v1(k)*rcb*cosinus-dx(k)*v1m)*dcosdcalpconst
+                dcosdcm(k) = ((dx(k)-v1(k))*v1m*rcb+ &
+                        cosinus*(dx(k)*v1m**2-v1(k)*rcb2))*dcosdcmconst
+             enddo
+             rcav = rcb/wcavsig
+             rcav11 = rcav**11
+             rcav12 = rcav11*rcav
+             constcav1 = 1-wcavchi*cos2
+             constcav2 = sqrt(constcav1)
+             constgb1 = 1/sqrt(1-wgbchi*cos2)
+             constgb2 = wgbeps*(1-wgbchip*cos2)**2
+             constdvan1 = 12*wvan1*wvan2**12*invrcs14
+             constdvan2 = 6*wvan1*wvan2**6*invrcs8
+!----------------------------------------------------------------------------
+!Gay-Berne term
+!---------------------------------------------------------------------------
+             sgb = 1/(1-constgb1+(rcb/wgbsig))
+             sgb6 = sgb**6
+             sgb7 = sgb6*sgb
+             sgb12 = sgb6**2
+             sgb13 = sgb12*sgb
+             Egb = constgb2*(sgb12-sgb6)
+             do k=1,3
+                dEgbdCat(k) = -constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+                 +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+     -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcat(k)
+                dEgbdCm(k) = constgb2/wgbsig*(12*sgb13-6*sgb7)*invrcb*dx(k) &
+                 +(constgb1**3*constgb2*wgbchi*cosinus*(12*sgb13-6*sgb7) &
+     -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcm(k)
+                dEgbdCalp(k) = (constgb1**3*constgb2*wgbchi*cosinus &
+                               *(12*sgb13-6*sgb7) &
+     -4*wgbeps*wgbchip*cosinus*(1-wgbchip*cos2)*(sgb12-sgb6))*dcosdcalp(k)
+             enddo
+!----------------------------------------------------------------------------
+!cavity term
+!---------------------------------------------------------------------------
+             cavnum = sqrt(rcav*constcav2)+wcav2*rcav*constcav2-wcav3
+             cavdenom = 1+wcav4*rcav12*constcav1**6
+             Ecav = wcav1*cavnum/cavdenom
+             invcavdenom2 = 1/cavdenom**2
+             dcavnumdcos = -wcavchi*cosinus/constcav2 &
+                    *(sqrt(rcav/constcav2)/2+wcav2*rcav)
+             dcavnumdr = (0.5*sqrt(constcav2/rcav)+wcav2*constcav2)/wcavsig
+             dcavdenomdcos = -12*wcav4*wcavchi*rcav12*constcav1**5*cosinus
+             dcavdenomdr = 12*wcav4/wcavsig*rcav11*constcav1**6
+             do k=1,3
+                dEcavdCat(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+     *dcosdcat(k)+(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+                dEcavdCm(k) = ((dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+     *dcosdcm(k)-(dcavnumdr*cavdenom-dcavdenomdr*cavnum)/rcb*dx(k))*wcav1*invcavdenom2
+                dEcavdCalp(k) = (dcavnumdcos*cavdenom-dcavdenomdcos*cavnum) &
+                             *dcosdcalp(k)*wcav1*invcavdenom2
+             enddo
+!----------------------------------------------------------------------------
+!van der Waals and dipole-charge interaction energy
+!---------------------------------------------------------------------------
+             Evan1 = wvan1*wvan2**12*invrcs12
+             do k=1,3
+                dEvan1Cat(k) = -v2(k)*constdvan1
+                dEvan1Cm(k) = 0.0d0
+                dEvan1Calp(k) = v2(k)*constdvan1
+             enddo
+             Evan2 = -wvan1*wvan2**6*invrcs6
+             do k=1,3
+                dEvan2Cat(k) = v2(k)*constdvan2
+                dEvan2Cm(k) = 0.0d0
+                dEvan2Calp(k) = -v2(k)*constdvan2
+             enddo
+             Edip = wdip1*cosinus*invrcb2-wdip2*(1-cos2)*invrcb4
+             do k=1,3
+                dEdipCat(k) = (-2*wdip1*cosinus*invrcb4 &
+                               +4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+                   +dcosdcat(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+                dEdipCm(k) = (2*wdip1*cosinus*invrcb4 &
+                             -4*wdip2*(1-cos2)*invrcb6)*dx(k) &
+                   +dcosdcm(k)*(wdip1*invrcb2+2*wdip2*cosinus*invrcb4)
+                dEdipCalp(k) = dcosdcalp(k)*(wdip1*invrcb2 &
+                                  +2*wdip2*cosinus*invrcb4)
+             enddo
+             if (energy_dec) write (iout,'(2i5,4(a6,f7.3))') i,j, &
+         ' E GB ',Egb,' ECav ',Ecav,' Evdw ',Evan1+Evan2,' Edip ',Edip
+             ecation_nucl=ecation_nucl+Ecav+Egb+Edip+Evan1+Evan2
+             do k=1,3
+                dEtotalCat(k) = dEcavdCat(k)+dEvan1Cat(k)+dEvan2Cat(k) &
+                                             +dEgbdCat(k)+dEdipCat(k)
+                dEtotalCm(k) = dEcavdCm(k)+dEvan1Cm(k)+dEvan2Cm(k) &
+                                           +dEgbdCm(k)+dEdipCm(k)
+                dEtotalCalp(k) = dEcavdCalp(k)+dEgbdCalp(k)+dEvan1Calp(k) &
+                                             +dEdipCalp(k)+dEvan2Calp(k)
+             enddo
+             do k=1,3
+                gg(k) = dEtotalCm(k)+dEtotalCalp(k)
+                gradnuclcatx(k,i)=gradnuclcatx(k,i)+dEtotalCm(k)
+                gradnuclcat(k,i)=gradnuclcat(k,i)+gg(k)
+                gradnuclcat(k,j)=gradnuclcat(k,j)+dEtotalCat(k)
+             enddo
+          enddo !j
+       enddo !i
+       return
+       end subroutine ecat_nucl
+
 !-----------------------------------------------------------------------------
 !-----------------------------------------------------------------------------
       subroutine eprot_sc_base(escbase)
@@ -23458,11 +23683,11 @@ chip1=chip(itypi)
          yj=c(2,j+nres)
          zj=c(3,j+nres)
       call to_box(xj,yj,zj)
-      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
@@ -24346,11 +24571,11 @@ chip1=chip(itypi)
          yj=(c(2,j)+c(2,j+1))/2.0
          zj=(c(3,j)+c(3,j+1))/2.0
      call to_box(xj,yj,zj)
-     call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
-      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
-      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
-       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!     call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+!      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
+!       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       xj=boxshift(xj-xi,boxxsize)
       yj=boxshift(yj-yi,boxysize)
       zj=boxshift(zj-zi,boxzsize)
@@ -25093,6 +25318,7 @@ chip1=chip(itypi)
          zj=c(3,j+nres)
      call to_box(xj,yj,zj)
      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+      write(iout,*) "KRUWA", i,j
       aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 &
@@ -26644,7 +26870,7 @@ chip1=chip(itypi)
 
 !c!-------------------------------------------------------------------
 !c! ecl
-       write(iout,*) "KURWA2",Rhead
+!       write(iout,*) "KURWA2",Rhead
        sparrow  = w1 * Qj * om1
        hawk     = w2 * Qj * Qj * (1.0d0 - sqom2)
        ECL = sparrow / Rhead**2.0d0 &
@@ -27345,7 +27571,13 @@ chip1=chip(itypi)
              yj=c(2,nres+j)
              zj=c(3,nres+j)
              call to_box(xj,yj,zj)
-             dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+!          call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+!          faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
+          dist_init=xj**2+yj**2+zj**2
+!             dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
 ! r_buff_list is a read value for a buffer 
              if (sqrt(dist_init).le.(r_cut_ele+r_buff_list)) then
 ! Here the list is created
@@ -27552,7 +27784,7 @@ chip1=chip(itypi)
 !      integer :: newcontlistppi(200*nres),newcontlistppj(200*nres)
       integer i,j,itypi,itypj,subchap,xshift,yshift,zshift,iint,ilist_pp,g_ilist_pp
       integer displ(0:nprocs),i_ilist_pp(0:nprocs),ierr
-            write(iout,*),"START make_pp",iatel_s,iatel_e,r_cut_ele+r_buff_list
+!            write(iout,*),"START make_pp",iatel_s,iatel_e,r_cut_ele+r_buff_list
             ilist_pp=0
       r_buff_list=5.0
       do i=iatel_s,iatel_e

diff --git a/source/unres/io.F90 b/source/unres/io.F90
index 86b40ef..e8b0579 100644 (file)
--- a/source/unres/io.F90
+++ b/source/unres/io.F90
@@ -807,6 +807,7 @@
       call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
       call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
       call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
+      call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
       call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
       call reada(weightcard,'WSCPHO',wscpho,0.0d0)
       call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
@@ -852,6 +853,7 @@
           weights(47)=wpepbase
           weights(48)=wscpho
           weights(49)=wpeppho
+          weights(50)=wcatnucl
 
       if(me.eq.king.or..not.out1file) &
        write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&

diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90
index dce2e69..61cb7b4 100644 (file)
--- a/source/unres/io_config.F90
+++ b/source/unres/io_config.F90
@@ -867,6 +867,9 @@
           enddo
         enddo
       endif
+
+
+
       if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
        if (oldion.eq.1) then
             do i=1,ntyp_molec(5)
@@ -882,6 +885,22 @@
             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 (iion,*) (vcatprm(k),k=1,ncatprotpram)
 !----------------------------------------------------
       allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
@@ -5337,6 +5356,8 @@
       open (itube,file=tubename,status='old')
       call getenv_loc('IONPAR',ionname)
       open (iion,file=ionname,status='old')
+      call getenv_loc('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old')
 #else
       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
         readonly)
@@ -5396,6 +5417,8 @@
       open (itube,file=tubename,status='old',action='read')
       call getenv_loc('IONPAR',ionname)
       open (iion,file=ionname,status='old',action='read')
+      call getenv_loc('IONPAR_NUCL',ionnuclname)
+      open (iionnucl,file=ionnuclname,status='old',action='read')
 
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)