wprowadzenie lipidow
authorAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 10 Feb 2015 11:54:18 +0000 (12:54 +0100)
committerAdam Sieradzan <adasko@piasek4.chem.univ.gda.pl>
Tue, 10 Feb 2015 11:54:18 +0000 (12:54 +0100)
14 files changed:
CMakeLists.txt
source/unres/src_MD-M/COMMON.CHAIN
source/unres/src_MD-M/COMMON.DERIV
source/unres/src_MD-M/COMMON.FFIELD
source/unres/src_MD-M/COMMON.INTERACT
source/unres/src_MD-M/COMMON.IOUNITS
source/unres/src_MD-M/COMMON.LOCAL
source/unres/src_MD-M/DIMENSIONS
source/unres/src_MD-M/energy_p_new-sep_barrier.F
source/unres/src_MD-M/energy_p_new_barrier.F
source/unres/src_MD-M/initialize_p.F
source/unres/src_MD-M/parmread.F
source/unres/src_MD-M/readrtns_CSA.F
source/unres/src_MD-M/unres.F

index 89bd8ef..bd9465a 100644 (file)
@@ -195,7 +195,7 @@ if(UNRES_NA_MMCE)
   add_subdirectory(source/cluster/unres/src)
   add_subdirectory(source/xdrfpdb/src)
   add_subdirectory(source/xdrfpdb/src-M)
-  add_subdirectory(source/maxlik/src_CSA)
+#  add_subdirectory(source/maxlik/src_CSA)
 else()
 
   add_subdirectory(source/unres/src_MD)
@@ -212,6 +212,6 @@ else()
   add_subdirectory(source/cluster/unres/src)
   add_subdirectory(source/xdrfpdb/src)
   add_subdirectory(source/xdrfpdb/src-M)
-  add_subdirectory(source/maxlik/src_CSA)
+#  add_subdirectory(source/maxlik/src_CSA)
 endif(UNRES_NA_MMCE)
 
index 777cc43..97cb82b 100644 (file)
@@ -14,5 +14,8 @@
      & nsup,nstart_sup,nstart_seq,
      & chain_length,iprzes,tabperm(maxperm,maxsym),nperm
       common /from_zscore/ nz_start,nz_end,iz_sc
-      double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad
+      double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
+     &buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
       common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad
+     &buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+
index 7688aeb..ba3f6c5 100644 (file)
@@ -1,5 +1,5 @@
-      double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gelc_long,
-     & gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp,
+      double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gelc_long
+     & gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp,gliptranc,gliptranx,
      & gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gloc_x,dtheta,dphi,dalpha,
      & domega,gscloc,gsclocx,gradcorr,gradcorr_long,gradcorr5_long,
      & gradcorr6_long,gcorr6_turn_long,gvdwx
@@ -7,21 +7,24 @@
       common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres),
      & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres),
      & gvdwc(3,maxres),gelc(3,maxres),gelc_long(3,maxres),
-     & gvdwpp(3,maxres),gvdwc_scpp(3,maxres),
+     & gvdwpp(3,maxres),gvdwc_scpp(3,maxres),gliptranc(3,maxres),
+     & gliptranx(3,maxres),
      & gradx_scp(3,maxres),gvdwc_scp(3,maxres),ghpbx(3,maxres),
-     & ghpbc(3,maxres),gloc(0:maxvar,2),gradcorr(3,maxres),
+     & ghpbc(3,maxres),gloc(maxvar,2),gradcorr(3,maxres),
      & gradcorr_long(3,maxres),gradcorr5_long(3,maxres),
      & gradcorr6_long(3,maxres),gcorr6_turn_long(3,maxres),
      & gradxorr(3,maxres),gradcorr5(3,maxres),gradcorr6(3,maxres),
      & gloc_x(maxvar,2),gel_loc(3,maxres),gel_loc_long(3,maxres),
      & gcorr3_turn(3,maxres),
-     & gcorr4_turn(3,maxres),gcorr6_turn(3,maxres),gradb(3,maxres),
+     &gcorr4_turn(3,maxres),gcorr6_turn(3,maxres),gradb(3,maxres),
      & gradbx(3,maxres),gel_loc_loc(maxvar),gel_loc_turn3(maxvar),
-     & gel_loc_turn4(maxvar),gel_loc_turn6(maxvar),gcorr_loc(maxvar),
+     & gel_loc_turn4(maxvar),gel_loc_turn6(maxvar),
+     & gcorr_loc(maxvar),
      & g_corr5_loc(maxvar),g_corr6_loc(maxvar),gsccorc(3,maxres),
      & gsccorx(3,maxres),gsccor_loc(maxres),dtheta(3,2,maxres),
      & gscloc(3,maxres),gsclocx(3,maxres),
-     & dphi(3,3,maxres),dalpha(3,3,maxres),domega(3,3,maxres),nfl,icg
+     & dphi(3,3,maxres),dalpha(3,3,maxres),domega(3,3,maxres),nfl,
+     & icg
       double precision derx,derx_turn
       common /deriv_loc/ derx(3,5,2),derx_turn(3,5,2)
       double precision dXX_C1tab(3,maxres),dYY_C1tab(3,maxres),
index d7d8cde..a63fe78 100644 (file)
@@ -6,7 +6,7 @@ C-----------------------------------------------------------------------
       integer n_ene_comp,rescale_mode
       common /ffield/ wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,
      &  wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
-     &  wturn6,wvdwpp,weights(n_ene),temp0,
+     &  wturn6,wvdwpp,weights(n_ene),wliptran,temp0,
      &  scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp,
      &  rescale_mode
       common /potentials/ potname(5)
index 37d3e88..4f7a354 100644 (file)
@@ -1,11 +1,14 @@
-      double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6
+      double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6,
+     &aa_lip,bb_lip,aa_aq,bb_aq
       integer expon,expon2
       integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,
      & ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr,iscpstart,
      & iscpend,iatsc_s,iatsc_e,
      & iatel_s,iatel_e,iatscp_s,iatscp_e,iatel_s_vdw,iatel_e_vdw,
      & ispp,iscp
-      common /interact/aa(ntyp,ntyp),bb(ntyp,ntyp),augm(ntyp,ntyp),
+      common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp),
+     & aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp),aa(ntyp,ntyp),bb(ntyp,ntyp),
+     & augm(ntyp,ntyp),
      & aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2),
      & expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr),
      & iend(maxres,maxint_gr),itype(maxres),itel(maxres),itypro,
@@ -31,3 +34,7 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters.
      & vbldsc0(maxbondterm,ntyp),akp,
      & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp),
      & distchainmax,nbondterm(ntyp)
+C 01/29/15 Lipidic parameters
+      double precision   pepliptran,liptranene
+      common /lipid/ pepliptran,liptranene(ntyp)
+
index 49b6db3..56e655e 100644 (file)
@@ -39,10 +39,10 @@ C CSA I/O units & files
 C Parameter files
       character*256 bondname,thetname,rotname,torname,tordname,
      &       fouriername,elename,sidename,scpname,sccorname,patname,
-     &       thetname_pdb,rotname_pdb
+     &       thetname_pdb,rotname_pdb,liptranname
       common /parfiles/ bondname,thetname,rotname,torname,tordname,
      &       fouriername,elename,sidename,scpname,sccorname,patname,
-     &       thetname_pdb,rotname_pdb
+     &       thetname_pdb,rotname_pdb,liptranname
       character*3 pot
 C-----------------------------------------------------------------------
 C INP    - main input file
index 9f0627b..5d1ced7 100644 (file)
@@ -40,7 +40,8 @@ C Virtual-bond lenghts
      & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end,
      & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end,
      & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,
-     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,
+     & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,ilip_start,
+     & ilip_end,
      & ibond_displ(0:max_fg_procs-1),ibond_count(0:max_fg_procs-1),
      & ithet_displ(0:max_fg_procs-1),ithet_count(0:max_fg_procs-1),
      & iphi_displ(0:max_fg_procs-1),iphi_count(0:max_fg_procs-1),
@@ -56,6 +57,6 @@ C Virtual-bond lenghts
      & iint_end,iphi1_start,iphi1_end,iint_count,iint_displ,ivec_displ,
      & ivec_count,iset_displ,itau_start,itau_end,
      & iset_count,ibond_displ,ibond_count,ithet_displ,ithet_count,
-     & iphi_displ,iphi_count,iphi1_displ,iphi1_count
+     & iphi_displ,iphi_count,iphi1_displ,iphi1_count,ilip_start,ilip_end
 C Inverses of the actual virtual bond lengths
       common /invlen/ vbld_inv(maxres2)
index 84cf3fd..fa5cf66 100644 (file)
@@ -6,17 +6,17 @@
 ********************************************************************************
 C Max. number of processors.
       integer maxprocs
-      parameter (maxprocs=2048)
+      parameter (maxprocs=1028)
 C Max. number of fine-grain processors
       integer max_fg_procs
 c      parameter (max_fg_procs=maxprocs)
-      parameter (max_fg_procs=512)
+      parameter (max_fg_procs=256)
 C Max. number of coarse-grain processors
       integer max_cg_procs
       parameter (max_cg_procs=maxprocs)
 C Max. number of AA residues
       integer maxres
-      parameter (maxres=800)
+      parameter (maxres=600)
 C Appr. max. number of interaction sites
       integer maxres2,maxres6,mmaxres2
       parameter (maxres2=2*maxres,maxres6=6*maxres)
@@ -95,7 +95,7 @@ C Max. number of conformations in the pool
       parameter (max_pool=10)
 C Number of energy components
       integer n_ene,n_ene2
-      parameter (n_ene=21,n_ene2=2*n_ene)
+      parameter (n_ene=22,n_ene2=2*n_ene)
 C Number of threads in deformation
       integer max_thread,max_thread2
       parameter (max_thread=4,max_thread2=2*max_thread)     
index 292bbee..bba2354 100644 (file)
@@ -6,7 +6,7 @@ C      if(r.lt.r_cut-rlamb) then
 C        sscale=1.0d0
 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
 C        gamm=(r-(r_cut-rlamb))/rlamb
-        sscale=1.0d0+r*r*(2*r-3.0d0)
+        sscalelip=1.0d0+r*r*(2*r-3.0d0)
 C      else
 C        sscale=0d0
 C      endif
@@ -20,7 +20,7 @@ C     if(r.lt.r_cut-rlamb) then
 C        sscagrad=0.0d0
 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
 C        gamm=(r-(r_cut-rlamb))/rlamb
-        sscagrad=r*(6*r-6.0d0)
+        sscagradlip=r*(6*r-6.0d0)
 C      else
 C        sscagrad=0.0d0
 C      endif
index a8c0323..7207b35 100644 (file)
@@ -199,6 +199,7 @@ c      print *,"Processor",myrank," computed UB"
 C
 C Calculate the SC local energy.
 C
+C      print *,"TU DOCHODZE?"
       call esc(escloc)
 c      print *,"Processor",myrank," computed USC"
 C
@@ -229,6 +230,7 @@ C
       else
         esccor=0.0d0
       endif
+C      print *,"PRZED MULIt"
 c      print *,"Processor",myrank," computed Usccorr"
 C 
 C 12/1/95 Multi-body terms
@@ -264,9 +266,11 @@ C  after the equilibration time
 C 01/27/2015 added by adasko
 C the energy component below is energy transfer into lipid environment 
 C based on partition function
+C      print *,"przed lipidami"
       if (wliptran.gt.0) then
-        call Eliptransfer
+        call Eliptransfer(eliptran)
       endif
+C      print *,"za lipidami"
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
 #endif
@@ -308,7 +312,7 @@ C
       energia(17)=estr
       energia(20)=Uconst+Uconst_back
       energia(21)=esccor
-      energia(22)=eliptrans
+      energia(22)=eliptran
 c    Here are the energies showed per procesor if the are more processors 
 c    per molecule then we sum it up in sum_energy subroutine 
 c      print *," Processor",myrank," calls SUM_ENERGY"
@@ -400,7 +404,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
-      energia(22)=eliptrans
+      eliptran=energia(22)
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
@@ -451,8 +455,8 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef MPI
       include 'mpif.h'
 #endif
-      double precision gradbufc(3,maxres),gradbufx(3,maxres),
-     &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
+      double precision gradbufc(3,0:maxres),gradbufx(3,0:maxres),
+     & glocbuf(4*maxres),gradbufc_sum(3,0:maxres),gloc_scbuf(3,0:maxres)
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
@@ -516,6 +520,8 @@ c      enddo
      &                wcorr6*gradcorr6_long(j,i)+
      &                wturn6*gcorr6_turn_long(j,i)+
      &                wstrain*ghpbc(j,i)
+     &                +wliptran*gliptranc(j,i)
+
         enddo
       enddo 
 #else
@@ -531,6 +537,7 @@ c      enddo
      &                wcorr6*gradcorr6_long(j,i)+
      &                wturn6*gcorr6_turn_long(j,i)+
      &                wstrain*ghpbc(j,i)
+     &                +wliptran*gliptranc(j,i)
         enddo
       enddo 
 #endif
@@ -982,6 +989,7 @@ C------------------------------------------------------------------------
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
+      eliptran=energia(22)
 #ifdef SPLITELE
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
      &  estr,wbond,ebe,wang,
@@ -990,7 +998,7 @@ C------------------------------------------------------------------------
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
      &  edihcnstr,ebr*nss,
-     &  Uconst,etot
+     &  Uconst,eliptran,wliptran,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -1014,6 +1022,7 @@ C------------------------------------------------------------------------
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
+     & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
 #else
       write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
@@ -1022,7 +1031,7 @@ C------------------------------------------------------------------------
      &  ecorr,wcorr,
      &  ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
      &  eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
-     &  ebr*nss,Uconst,etot
+     &  ebr*nss,Uconst,eliptran,wliptran,etot
    10 format (/'Virtual-chain energies:'//
      & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
      & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
@@ -1045,6 +1054,7 @@ C------------------------------------------------------------------------
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST=',1pE16.6,' (Constraint energy)'/ 
+     & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
      & 'ETOT=  ',1pE16.6,' (total)')
 #endif
       return
@@ -1099,6 +1109,7 @@ C Change 12/1/95 to calculate four-body interactions
 c           write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
             eps0ij=eps(itypi,itypj)
             fac=rrij**expon2
+C have you changed here?
             e1=fac*fac*aa(itypi,itypj)
             e2=fac*bb(itypi,itypj)
             evdwij=e1+e2
@@ -1249,6 +1260,7 @@ C
             rij=1.0D0/r_inv_ij 
             r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
             fac=r_shift_inv**expon
+C have you changed here?
             e1=fac*fac*aa(itypi,itypj)
             e2=fac*bb(itypi,itypj)
             evdwij=e_augm+e1+e2
@@ -1376,6 +1388,7 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives.
             call sc_angular
 C Calculate whole angle-dependent part of epsilon and contributions
 C to its derivatives
+C have you changed here?
             fac=(rrij*sigsq)**expon2
             e1=fac*fac*aa(itypi,itypj)
             e2=fac*bb(itypi,itypj)
@@ -1636,6 +1649,8 @@ cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
             rij_shift=1.0D0/rij_shift 
             fac=rij_shift**expon
+C here to start with
+C            if (c(i,3).gt.
             e1=fac*fac*aa(itypi,itypj)
             e2=fac*bb(itypi,itypj)
             evdwij=eps1*eps2rt*eps3rt*(e1+e2)
@@ -1674,6 +1689,7 @@ C Calculate the radial part of the gradient
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
             call sc_grad
+            endif
             ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
@@ -2982,6 +2998,8 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
 C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
       do i=iturn3_start,iturn3_end
+        if (i.le.1) cycle
+C        write(iout,*) "tu jest i",i
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i+3).eq.ntyp1
@@ -3009,6 +3027,7 @@ C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
         num_cont_hb(i)=num_conti
       enddo
       do i=iturn4_start,iturn4_end
+        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
@@ -3071,6 +3090,7 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
       do i=iatel_s,iatel_e
+        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i-1).eq.ntyp1
@@ -3123,7 +3143,8 @@ c        endif
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
-c          write (iout,*) i,j,itype(i),itype(j)
+C          write (iout,*) i,j
+         if (j.le.1) cycle
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
      & .or.itype(j+2).eq.ntyp1
      & .or.itype(j-1).eq.ntyp1
@@ -4660,15 +4681,14 @@ c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
 c     &    dhpb(i),dhpb1(i),forcon(i)
 C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
-        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
-     & iabs(itype(jjj)).eq.1) then
+C        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C     & iabs(itype(jjj)).eq.1) then
 cmc        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
 C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
         if (.not.dyn_ss .and. i.le.nss) then
 C 15/02/13 CC dynamic SSbond - additional check
          if (ii.gt.nres 
      &       .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then 
->>>>>>> f5379d3246c4bd95e946c4d35d4a1c13e329c4cb
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
          endif
@@ -9674,6 +9694,7 @@ crc      print *,((prod(i,j),i=1,2),j=1,2)
       end
 CCC----------------------------------------------
       subroutine Eliptransfer(eliptran)
+      implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
       include 'COMMON.GEO'
       include 'COMMON.VAR'
@@ -9687,6 +9708,7 @@ CCC----------------------------------------------
       include 'COMMON.CONTROL'
       include 'COMMON.SPLITELE'
       include 'COMMON.SBRIDGE'
+C      print *,"wchodze"
 C structure of box:
 C      water
 C--bordliptop-- buffore starts
@@ -9695,67 +9717,91 @@ C      lipid
 C--buflipbot--- lipid ends buffore starts
 C--bordlipbot--buffore ends
       eliptran=0.0
-      do i=1,nres
+      do i=ilip_start,ilip_end
+        if (itype(i).eq.ntyp1) cycle
+
+        positi=(mod((c(3,i)+c(3,i+1)),boxzsize))
+        if (positi.le.0) positi=positi+boxzsize
+C        print *,i
 C first for peptide groups
 c for each residue check if it is in lipid or lipid water border area
-       if ((mod(c(3,i),boxzsize).gt.bordlipbot)
-     &.and.(mod(c(3,i),boxzsize).lt.bordliptop)) then
+       if ((positi.gt.bordlipbot)
+     &.and.(positi.lt.bordliptop)) then
 C the energy transfer exist
-        if (mod(c(3,i),boxzsize).lt.buflipbot) then
+        if (positi.lt.buflipbot) then
 C what fraction I am in
          fracinbuf=1.0d0-
-     &        ((mod(c(3,i),boxzsize)-bordlipbot)/lipbufthick)
+     &        ((positi-bordlipbot)/lipbufthick)
 C lipbufthick is thickenes of lipid buffore
-         ssslip=sscale(fracinbuf)
+         sslip=sscalelip(fracinbuf)
          ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
-         eliptran=eliptran+sslip
-         gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+         eliptran=eliptran+sslip*pepliptran
+         gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0
+         gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0
+C         gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
+
 C         print *,"doing sccale for lower part"
-        elseif (mod(c(3,i),boxzsize).gt.bufliptop) then
-         fracinbuf=1.0d0-((bordliptop-mod(c(3,i),boxzsize))/lipbufthick)
-         ssslip=sscale(fracinbuf)
+        elseif (positi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+         sslip=sscalelip(fracinbuf)
          ssgradlip=sscagradlip(fracinbuf)/lipbufthick
-         eliptran=eliptran+sslip
-         gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+         eliptran=eliptran+sslip*pepliptran
+         gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+         gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C         gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
           print *, "doing sscalefor top part"
         else
-         eliptran=eliptran+1.0d0
+         eliptran=eliptran+pepliptran
          print *,"I am in true lipid"
         endif
 C       else
 C       eliptran=elpitran+0.0 ! I am in water
        endif
        enddo
+C       print *, "nic nie bylo w lipidzie?"
 C now multiply all by the peptide group transfer factor
-       eliptran=eliptran*pepliptran
+C       eliptran=eliptran*pepliptran
 C now the same for side chains
-       do i=1,nres
+       do i=ilip_start,ilip_end
+        if (itype(i).eq.ntyp1) cycle
+        positi=(mod(c(3,i+nres),boxzsize))
+        if (positi.le.0) positi=positi+boxzsize
+C       print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
 c for each residue check if it is in lipid or lipid water border area
-       if ((mod(c(3,i+nres),boxzsize).gt.bordlipbot)
-     & .and.(mod(c(3,i+nres),boxzsize).lt.bordliptop)) then
+C       respos=mod(c(3,i+nres),boxzsize)
+       if ((positi.gt.bordlipbot)
+     & .and.(positi.lt.bordliptop)) then
 C the energy transfer exist
-        if (mod(c(3,i+nres),boxzsize).lt.buflipbot) then
+        if (positi.lt.buflipbot) then
          fracinbuf=1.0d0-
-     &     ((mod(c(3,i+nres),boxzsize)-bordlipbot)/lipbufthick)
+     &     ((positi-bordlipbot)/lipbufthick)
 C lipbufthick is thickenes of lipid buffore
-         ssslip=sscale(fracinbuf)
+         sslip=sscalelip(fracinbuf)
          ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
          eliptran=eliptran+sslip*liptranene(itype(i))
-         gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
+         gliptranx(3,i)=gliptranx(3,i)
+     &+ssgradlip*liptranene(itype(i))/2.0d0
+         gliptranc(3,i-1)=
+     &+ssgradlip*liptranene(itype(i))
          print *,"doing sccale for lower part"
-        elseif (mod(c(3,i+nres),boxzsize).gt.bufliptop) then
+        elseif (positi.gt.bufliptop) then
          fracinbuf=1.0d0-
-     &((bordliptop-mod(c(3,i+nres),boxzsize))/lipbufthick)
-         ssslip=sscale(fracinbuf)
+     &((bordliptop-positi)/lipbufthick)
+         sslip=sscalelip(fracinbuf)
          ssgradlip=sscagradlip(fracinbuf)/lipbufthick
          eliptran=eliptran+sslip*liptranene(itype(i))
-         gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
-          print *, "doing sscalefor top part"
+         gliptranx(3,i)=gliptranx(3,i)
+     &+ssgradlip*liptranene(itype(i))/2.0d0
+         gliptranc(3,i-1)=
+     &+ssgradlip*liptranene(itype(i))
+          print *, "doing sscalefor top part",sslip,fracinbuf
         else
          eliptran=eliptran+liptranene(itype(i))
          print *,"I am in true lipid"
         endif
+        endif ! if in lipid or buffor
 C       else
 C       eliptran=elpitran+0.0 ! I am in water
        enddo
-
+       return
+       end
index 818abde..c4baab2 100644 (file)
@@ -119,6 +119,14 @@ C
       icsa_in=40
 crc for ifc error 118
       icsa_pdb=42
+C Lipidic input file for parameters range 60-79
+      iliptranpar=60
+C input file for transfer sidechain and peptide group inside the
+C lipidic environment if lipid is implicite
+
+C DNA input files for parameters range 80-99
+C Suger input files for parameters range 100-119
+C All-atom input files for parameters range 120-149
 C
 C Set default weights of the energy terms.
 C
@@ -632,6 +640,8 @@ C Partition local interactions
       call int_bounds(nct-nnt,ibondp_start,ibondp_end) 
       ibondp_start=ibondp_start+nnt
       ibondp_end=ibondp_end+nnt
+      call int_bounds(nres,ilip_start,ilip_end)
+      ilip_start=ilip_start
       call int_bounds1(nres-1,ivec_start,ivec_end) 
 c      print *,"Processor",myrank,fg_rank,fg_rank1,
 c     &  " ivec_start",ivec_start," ivec_end",ivec_end
@@ -1156,6 +1166,8 @@ c        write (iout,*) "MPI_ROTAT2",MPI_ROTAT2
       iset_end=nres+1
       iint_start=2
       iint_end=nres-1
+      ilip_start=1
+      ilip_end=nres
 #endif
       return
       end 
index bd2165b..0fc1ffa 100644 (file)
@@ -97,6 +97,12 @@ c
           enddo
         enddo
       endif
+C reading lipid parameters
+       read(iliptranpar,*) pepliptran
+       do i=1,ntyp
+       read(iliptranpar,*) liptranene(i)
+       enddo
+       close(iliptranpar)
 #ifdef CRYST_THETA
 C
 C Read the parameters of the probability distribution/energy expression 
index 68d63c4..9c34d47 100644 (file)
@@ -224,6 +224,19 @@ C Reading the dimensions of box in x,y,z coordinates
 c Cutoff range for interactions
       call reada(controlcard,"R_CUT",r_cut,15.0d0)
       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+      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
+C      endif
+      if ((bordliptop.gt.boxzsize).or.(borlipbot.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
       if (me.eq.king .or. .not.out1file ) 
      &  write (iout,*) "DISTCHAINMAX",distchainmax
       
@@ -586,6 +599,7 @@ C Read weights of the subsequent energy terms.
        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
        call reada(weightcard,'TEMP0',temp0,300.0d0)
+       call reada(weightcard,'WLT',wliptran,0.0D0)
        if (index(weightcard,'SOFT').gt.0) ipot=6
 C 12/1/95 Added weight for the multi-body term WCORR
        call reada(weightcard,'WCORRH',wcorr,1.0D0)
@@ -1972,6 +1986,8 @@ C Get parameter filenames and open the parameter files.
       open (ielep,file=elename,status='old',readonly)
       call getenv_loc('SIDEPAR',sidename)
       open (isidep,file=sidename,status='old',readonly)
+      call getenv_loc('LIPTRANPAR',liptranname)
+      open (iliptranpar,file=liptranname,status='old',action='read')
 #ifndef CRYST_SC
       call getenv_loc('ROTPARPDB',rotname_pdb)
       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
index 4903ec0..642b361 100644 (file)
@@ -196,23 +196,12 @@ c---------------------------------------------------------------------------
       double precision energy(0:n_ene)
       double precision energy_long(0:n_ene),energy_short(0:n_ene)
       double precision varia(maxvar)
-<<<<<<< HEAD
-      if (indpdb.eq.0)     call chainbuild
-      time00=MPI_Wtime()
-      print *,'dc',c(1,1)
-      if (indpdb.ne.0) then
-      dc(1,0)=c(1,1)
-      dc(2,0)=c(2,1)
-      dc(3,0)=c(3,1)
-      endif
-=======
       if (indpdb.eq.0) call chainbuild
 #ifdef MPI
       time00=MPI_Wtime()
 #else
       time00=tcpu()
 #endif
->>>>>>> f5379d3246c4bd95e946c4d35d4a1c13e329c4cb
       call chainbuild_cart
       print *,'dc',dc(1,0),dc(2,0),dc(3,0)
       if (split_ene) then