From 0c5fb524ca994d617a32adfb678302e2f3dcef1e Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Tue, 10 Feb 2015 12:54:18 +0100 Subject: [PATCH] wprowadzenie lipidow --- CMakeLists.txt | 4 +- source/unres/src_MD-M/COMMON.CHAIN | 5 +- source/unres/src_MD-M/COMMON.DERIV | 17 +-- source/unres/src_MD-M/COMMON.FFIELD | 2 +- source/unres/src_MD-M/COMMON.INTERACT | 11 +- source/unres/src_MD-M/COMMON.IOUNITS | 4 +- source/unres/src_MD-M/COMMON.LOCAL | 5 +- source/unres/src_MD-M/DIMENSIONS | 8 +- source/unres/src_MD-M/energy_p_new-sep_barrier.F | 4 +- source/unres/src_MD-M/energy_p_new_barrier.F | 124 +++++++++++++++------- source/unres/src_MD-M/initialize_p.F | 12 +++ source/unres/src_MD-M/parmread.F | 6 ++ source/unres/src_MD-M/readrtns_CSA.F | 16 +++ source/unres/src_MD-M/unres.F | 11 -- 14 files changed, 156 insertions(+), 73 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 89bd8ef..bd9465a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 777cc43..97cb82b 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -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 + diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 7688aeb..ba3f6c5 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -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), diff --git a/source/unres/src_MD-M/COMMON.FFIELD b/source/unres/src_MD-M/COMMON.FFIELD index d7d8cde..a63fe78 100644 --- a/source/unres/src_MD-M/COMMON.FFIELD +++ b/source/unres/src_MD-M/COMMON.FFIELD @@ -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) diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT index 37d3e88..4f7a354 100644 --- a/source/unres/src_MD-M/COMMON.INTERACT +++ b/source/unres/src_MD-M/COMMON.INTERACT @@ -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) + diff --git a/source/unres/src_MD-M/COMMON.IOUNITS b/source/unres/src_MD-M/COMMON.IOUNITS index 49b6db3..56e655e 100644 --- a/source/unres/src_MD-M/COMMON.IOUNITS +++ b/source/unres/src_MD-M/COMMON.IOUNITS @@ -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 diff --git a/source/unres/src_MD-M/COMMON.LOCAL b/source/unres/src_MD-M/COMMON.LOCAL index 9f0627b..5d1ced7 100644 --- a/source/unres/src_MD-M/COMMON.LOCAL +++ b/source/unres/src_MD-M/COMMON.LOCAL @@ -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) diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index 84cf3fd..fa5cf66 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -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) diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index 292bbee..bba2354 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -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 diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index a8c0323..7207b35 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -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 diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 818abde..c4baab2 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -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 diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index bd2165b..0fc1ffa 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -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 diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 68d63c4..9c34d47 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -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') diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F index 4903ec0..642b361 100644 --- a/source/unres/src_MD-M/unres.F +++ b/source/unres/src_MD-M/unres.F @@ -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 -- 1.7.9.5