From affbf55b871a9bcc8fdb91efeb8e6bf7a3d8d003 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Fri, 14 Jul 2017 08:57:59 +0200 Subject: [PATCH] working nano --- PARAM/TiO2_fin.parm | 25 ++ source/unres/control.F90 | 6 + source/unres/data/control_data.f90 | 3 +- source/unres/data/energy_data.f90 | 9 +- source/unres/data/io_units.f90 | 2 +- source/unres/data/names.f90 | 11 +- source/unres/energy.f90 | 685 +++++++++++++++++++++++++++++++++++- source/unres/io.f90 | 1 + source/unres/io_config.f90 | 56 ++- 9 files changed, 773 insertions(+), 25 deletions(-) create mode 100644 PARAM/TiO2_fin.parm diff --git a/PARAM/TiO2_fin.parm b/PARAM/TiO2_fin.parm new file mode 100644 index 0000000..a03bc3b --- /dev/null +++ b/PARAM/TiO2_fin.parm @@ -0,0 +1,25 @@ + 3.929 3.894 -241.329 50.844 288.017 6.165E-10 0.0 ! Pep + 0.174 4.496 2.222 -0.393 -3.104 -1.392E-12 0.0 ! Cys + 4.173 4.136 -143.035 28.494 180.705 1.086E-10 0.0 ! Met + 10.793 4.152 -250.420 48.235 328.091 7.237E-11 0.0 ! Phe + 9.872 3.799 -308.062 62.468 382.714 3.262E-10 0.0 ! Ile + 12.025 4.056 -448.084 90.596 559.212 2.482E-10 0.0 ! Leu + 3.929 3.894 -241.329 50.844 288.017 6.165E-10 0.0 ! Val + 0.285 5.797 738.577 -155.331 -884.672 3.662E-09 0.0 ! Trp + 12.085 4.158 -275.901 53.658 358.248 1.287E-10 0.0 ! Tyr + 0.142 4.512 2.390 -0.478 -2.867 -3.237E-12 0.0 ! Ala + 0.000 0.000 0.000 0.000 0.000 0.000E+00 0.0 ! Gly + 4.779 2.663 -43.436 9.048 51.955 1.735E-10 0.0 ! Thr + 10.487 2.236 -42.090 7.652 56.064 2.307E-11 0.0 ! Ser + 1.582 3.601 673.244 -173.206 -646.238 4.709E-08 0.0 ! Gln + 9.692 3.299 -610.854 141.974 665.900 1.676E-08 0.0 ! Asn + 5.499 3.891 -353.183 80.947 391.185 4.071E-09 0.0 ! Glu + 9.527 2.826 -485.295 112.978 523.938 2.594E-08 0.0 ! Asp + 6.267 3.952 -176.750 35.322 222.723 1.641E-10 0.0 ! His + 0.755 4.649 6.176 -0.945 -9.823 1.249E-13 0.0 ! Arg + 0.112 5.414 8.877 -1.471 -13.249 -2.915E-13 0.0 ! Lys + 4.954 4.214 -188.800 37.556 238.429 1.111E-10 0.0 ! Pro + 0.000 0.000 0.000 0.000 0.000 0.000E+00 0.00 ! Pro + 0.000 0.000 0.000 0.000 0.000 0.000E+00 0.00 ! Pro + 0.000 0.000 0.000 0.000 0.000 0.000E+00 0.00 ! Pro + 0.000 0.000 0.000 0.000 0.000 0.000E+00 0.00 ! Pro diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 6a8befd..f45a88e 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -230,6 +230,7 @@ izsc=32 #endif iliptranpar=60 + itube=61 #if defined(WHAM_RUN) || defined(CLUSTER) ! ! setting the mpi variables for WHAM @@ -799,6 +800,9 @@ call int_bounds(nres,ilip_start,ilip_end) ilip_start=ilip_start ilip_end=ilip_end + call int_bounds(nres-1,itube_start,itube_end) + itube_start=itube_start + itube_end=itube_end if (ndih_constr.eq.0) then idihconstr_start=1 idihconstr_end=0 @@ -1326,6 +1330,8 @@ iint_end=nres-1 ilip_start=1 ilip_end=nres + itube_start=1 + itube_end=nres #endif !el common /przechowalnia/ ! deallocate(iturn3_start_all) diff --git a/source/unres/data/control_data.f90 b/source/unres/data/control_data.f90 index ddad045..1b4e429 100644 --- a/source/unres/data/control_data.f90 +++ b/source/unres/data/control_data.f90 @@ -29,7 +29,8 @@ ! commom.control ! common /cntrl/ integer :: modecalc,iscode,indpdb,indback,indphi,iranconf,& - icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr,shield_mode + icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr,shield_mode,& + tubemode logical :: minim,refstr,pdbref,overlapsc,& energy_dec,sideadd,lsecondary,read_cart,unres_pdb,& vdisulf,searchsc,lmuca,dccart,extconf,out1file,& diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index 65bf013..295d42a 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -64,7 +64,7 @@ integer :: rescale_mode real(kind=8) :: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang,& wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,& - wturn6,wvdwpp,wliptran,wshield,lipscale + wturn6,wvdwpp,wliptran,wshield,lipscale,wtube #ifdef CLUSTER real(kind=8) :: scalscp #endif @@ -89,6 +89,11 @@ ! common.interact ! common /interact/ real(kind=8),dimension(:,:),allocatable :: aa_aq,bb_aq,augm,aa_lip,bb_lip !(ntyp,ntyp) + real(kind=8),dimension(:),allocatable :: sc_aa_tube_par,sc_bb_tube_par,& + acavtub,bcavtub,ccavtub,dcavtub,tubetranene + real(kind=8) :: acavtubpep,bcavtubpep,ccavtubpep,dcavtubpep, & + tubetranenepep,pep_aa_tube,pep_bb_tube,tubeR0 + real(kind=8),dimension(3) :: tubecenter real(kind=8),dimension(:,:),allocatable :: aad,bad !(ntyp,2) real(kind=8),dimension(2,2) :: app,bpp,ael6,ael3 integer :: expon,expon2, nnt,nct,itypro @@ -150,7 +155,7 @@ 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,& - ilip_start,ilip_end + ilip_start,ilip_end,itube_start,itube_end integer,dimension(:),allocatable :: ibond_displ,ibond_count,& ithet_displ,ithet_count,iphi_displ,iphi_count,iphi1_displ,& iphi1_count,ivec_displ,ivec_count,iset_displ,iset_count,& diff --git a/source/unres/data/io_units.f90 b/source/unres/data/io_units.f90 index 1a88d8c..ebda347 100644 --- a/source/unres/data/io_units.f90 +++ b/source/unres/data/io_units.f90 @@ -43,7 +43,7 @@ ! common /parfiles/ character(len=256) :: bondname,thetname,rotname,torname,tordname,& fouriername,elename,sidename,scpname,sccorname,patname,& - thetname_pdb,rotname_pdb,liptranname + thetname_pdb,rotname_pdb,liptranname,tubename !----------------------------------------------------------------------- ! INP - main input file ! IOUT - list file diff --git a/source/unres/data/names.f90 b/source/unres/data/names.f90 index d5a23a2..7885f0b 100644 --- a/source/unres/data/names.f90 +++ b/source/unres/data/names.f90 @@ -28,7 +28,7 @@ !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Number of energy components - integer,parameter :: n_ene=21 + integer,parameter :: n_ene=25 integer :: n_ene2=2*n_ene !----------------------------------------------------------------------------- ! common.names @@ -51,15 +51,18 @@ (/"EVDW SC-SC","EVDW2 SC-p","EES p-p ","ECORR4 ","ECORR5 ",& "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",& "EBE bend ","ESC SCloc ","ETORS ","ETORSD ","EHPB ","EVDWPP ",& - "EVDW2_14 ","ESTR ","ESCCOR ","EDIHC ","EVDW_T "/) + "EVDW2_14 ","ESTR ","ESCCOR ","EDIHC ","EVDW_T ",& + "ELT "," "," ","ETUBE " /) + character(len=10),dimension(n_ene) :: wname = & (/"WSC ","WSCP ","WELEC" ,"WCORR ","WCORR5 ","WCORR6 ","WEL_LOC ",& "WTURN3 ","WTURN4 ","WTURN6 ","WANG ","WSCLOC ","WTOR ","WTORD ",& - "WHPB ","WVDWPP ","WSCP14 ","WBOND ","WSCCOR ","WDIHC ","WSC "/) + "WHPB ","WVDWPP ","WSCP14 ","WBOND ","WSCCOR ","WDIHC ","WSC ",& + "WLT "," "," ","WTUBE " /) 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/) + (/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/) !#endif !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 77c46dd..c166edf 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -123,7 +123,7 @@ gshieldc_ec,gshieldc_loc_ec,gshieldx_t3, & gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, & gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,& - grad_shield !(3,maxres) + grad_shield,gg_tube,gg_tube_sc !(3,maxres) ! real(kind=8),dimension(:,:),allocatable :: gloc,gloc_x !(maxvar,2) real(kind=8),dimension(:,:),allocatable :: gel_loc,gel_loc_long,& gcorr3_turn,gcorr4_turn,gcorr6_turn,gradb,gradbx !(3,maxres) @@ -222,7 +222,7 @@ integer :: n_corr,n_corr1,ierror real(kind=8) :: etors,edihcnstr,etors_d,esccor,ehpb real(kind=8) :: evdw,evdw1,evdw2,evdw2_14,escloc,ees,eel_loc - real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran + real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6 #ifdef MPI @@ -520,6 +520,15 @@ else eliptran=0.0d0 endif + if (tubemode.eq.1) then + call calctube(etube) + else if (tubemode.eq.2) then + call calctube2(etube) + elseif (tubemode.eq.3) then + call calcnano(etube) + else + etube=0.0d0 + endif #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 @@ -563,6 +572,7 @@ energia(20)=Uconst+Uconst_back energia(21)=esccor energia(22)=eliptran + energia(25)=etube ! Here are the energies showed per procesor if the are more processors ! per molecule then we sum it up in sum_energy subroutine ! print *," Processor",myrank," calls SUM_ENERGY" @@ -604,7 +614,7 @@ real(kind=8) :: evdw,evdw2,evdw2_14,ees,evdw1,ecorr,ecorr5,ecorr6 real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot, & - eliptran + eliptran,etube integer :: i #ifdef MPI integer :: ierr @@ -665,20 +675,21 @@ Uconst=energia(20) esccor=energia(21) eliptran=energia(22) + etube=energia(25) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & - +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & - +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+wtube*etube #endif energia(0)=etot ! detecting NaNQ @@ -801,7 +812,8 @@ !el local variables real(kind=8) :: etot,evdw,evdw2,ees,evdw1,ecorr,ecorr5,ecorr6,eel_loc real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc - real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran + real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,& + etube etot=energia(0) evdw=energia(1) @@ -832,7 +844,7 @@ Uconst=energia(20) esccor=energia(21) eliptran=energia(22) - + etube=energia(25) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -841,7 +853,7 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,& edihcnstr,ebr*nss,& - Uconst,eliptran,wliptran,etot + Uconst,eliptran,wliptran,etube,wtube,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -866,6 +878,7 @@ 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/& + 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -874,7 +887,7 @@ ecorr,wcorr,& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,& - ebr*nss,Uconst,eliptran,wliptran,etot + ebr*nss,Uconst,eliptran,wliptran,etube,wtube,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -898,6 +911,7 @@ 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & + 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -10131,7 +10145,9 @@ +wcorr*gshieldc_ec(j,i) & +wturn3*gshieldc_t3(j,i)& +wturn4*gshieldc_t4(j,i)& - +wel_loc*gshieldc_ll(j,i) + +wel_loc*gshieldc_ll(j,i)& + +wtube*gg_tube(j,i) + enddo @@ -10153,7 +10169,9 @@ +welec*gshieldc(j,i)& +wcorr*gshieldc_ec(j,i) & +wturn4*gshieldc_t4(j,i) & - +wel_loc*gshieldc_ll(j,i) + +wel_loc*gshieldc_ll(j,i)& + +wtube*gg_tube(j,i) + enddo @@ -10309,7 +10327,9 @@ +wturn4*gshieldc_t4(j,i) & +wturn4*gshieldc_loc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) & - +wel_loc*gshieldc_loc_ll(j,i) + +wel_loc*gshieldc_loc_ll(j,i) & + +wtube*gg_tube(j,i) + #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & @@ -10340,7 +10360,9 @@ +wturn4*gshieldc_t4(j,i) & +wturn4*gshieldc_loc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) & - +wel_loc*gshieldc_loc_ll(j,i) + +wel_loc*gshieldc_loc_ll(j,i) & + +wtube*gg_tube(j,i) + #endif @@ -10354,7 +10376,9 @@ +wcorr*gshieldx_ec(j,i) & +wturn3*gshieldx_t3(j,i) & +wturn4*gshieldx_t4(j,i) & - +wel_loc*gshieldx_ll(j,i) + +wel_loc*gshieldx_ll(j,i)& + +wtube*gg_tube_sc(j,i) + enddo enddo @@ -15860,7 +15884,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gshieldx_ll(j,i)=0.0d0 gshieldc_ll(j,i)=0.0d0 gshieldc_loc_ll(j,i)=0.0d0 - + gg_tube(j,i)=0.0d0 + gg_tube_sc(j,i)=0.0d0 do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo @@ -17904,6 +17929,634 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' enddo return end subroutine Eliptransfer +!----------------------------------NANO FUNCTIONS +!C----------------------------------------------------------------------- +!C----------------------------------------------------------- +!C This subroutine is to mimic the histone like structure but as well can be +!C utilizet to nanostructures (infinit) small modification has to be used to +!C make it finite (z gradient at the ends has to be changes as well as the x,y +!C gradient has to be modified at the ends +!C The energy function is Kihara potential +!C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +!C 4eps is depth of well sigma is r_minimum r is distance from center of tube +!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +!C simple Kihara potential + subroutine calctube(Etube) + real(kind=8) :: vectube(3),enetube(nres*2) + real(kind=8) :: Etube,xtemp,xminact,yminact,& + ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi, & + sc_aa_tube,sc_bb_tube + integer :: i,j,iti + Etube=0.0d0 + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=0.0d0 + enddo +!C first we calculate the distance from tube center +!C for UNRES + do i=itube_start,itube_end +!C lets ommit dummy atoms for now + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +!C now calculate distance from center of tube and direction vectors + xmin=boxxsize + ymin=boxysize +! Find minimum distance in periodic box + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) + +!C as the tube is infinity we do not calculate the Z-vector use of Z +!C as chosen axis + vectube(3)=0.0d0 +!C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +!C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +!C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +!C and its 6 power + rdiff6=rdiff**6.0d0 +!C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6 +!C write(iout,*) "TU13",i,rdiff6,enetube(i) +!C print *,rdiff,rdiff6,pep_aa_tube +!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +!C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6- & + 6.0d0*pep_bb_tube)/rdiff6/rdiff +!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +!C &rdiff,fac +!C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + enddo +!C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) +!C print *,gg_tube(1,0),"TU" + + + do i=itube_start,itube_end +!C Lets not jump over memory as we use many times iti + iti=itype(i) +!C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) & +!C in UNRES uncomment the line below as GLY has no side-chain... +!C .or.(iti.eq.10) + ) cycle + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp +!C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +!C & tubecenter(2) + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +!C as the tube is infinity we do not calculate the Z-vector use of Z +!C as chosen axis + vectube(3)=0.0d0 +!C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +!C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r + +!C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +!C and its 6 power + rdiff6=rdiff**6.0d0 +!C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 + fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- & + 6.0d0*sc_bb_tube/rdiff6/rdiff +!C now direction of gg_tube vector + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + enddo + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + enddo +!C print *,"ETUBE", etube + return + end subroutine calctube +!C TO DO 1) add to total energy +!C 2) add to gradient summation +!C 3) add reading parameters (AND of course oppening of PARAM file) +!C 4) add reading the center of tube +!C 5) add COMMONs +!C 6) add to zerograd +!C 7) allocate matrices + + +!C----------------------------------------------------------------------- +!C----------------------------------------------------------- +!C This subroutine is to mimic the histone like structure but as well can be +!C utilizet to nanostructures (infinit) small modification has to be used to +!C make it finite (z gradient at the ends has to be changes as well as the x,y +!C gradient has to be modified at the ends +!C The energy function is Kihara potential +!C E=4esp*((sigma/(r-r0))^12 - (sigma/(r-r0))^6) +!C 4eps is depth of well sigma is r_minimum r is distance from center of tube +!C and r0 is the excluded size of nanotube (can be set to 0 if we want just a +!C simple Kihara potential + subroutine calctube2(Etube) + real(kind=8) :: vectube(3),enetube(nres*2) + real(kind=8) :: Etube,xtemp,xminact,yminact,& + ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,positi,fracinbuf,& + sstube,ssgradtube,sc_aa_tube,sc_bb_tube + integer:: i,j,iti + Etube=0.0d0 + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=0.0d0 + enddo +!C first we calculate the distance from tube center +!C first sugare-phosphate group for NARES this would be peptide group +!C for UNRES + do i=itube_start,itube_end +!C lets ommit dummy atoms for now + + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +!C now calculate distance from center of tube and direction vectors +!C vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) +!C if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize +!C vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) +!C if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + +!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) + +!C as the tube is infinity we do not calculate the Z-vector use of Z +!C as chosen axis + vectube(3)=0.0d0 +!C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +!C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +!C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +!C and its 6 power + rdiff6=rdiff**6.0d0 +!C THIS FRAGMENT MAKES TUBE FINITE + positi=mod((c(3,i)+c(3,i+1))/2.0d0,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 +!C respos=mod(c(3,i+nres),boxzsize) +!C print *,positi,bordtubebot,buftubebot,bordtubetop + if ((positi.gt.bordtubebot) & + .and.(positi.lt.bordtubetop)) then +!C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- & + ((positi-bordtubebot)/tubebufthick) +!C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick +!C print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i)=enetube(i)+sstube*tubetranenepep +!C gg_tube_SC(3,i)=gg_tube_SC(3,i) +!C &+ssgradtube*tubetranene(itype(i)) +!C gg_tube(3,i-1)= gg_tube(3,i-1) +!C &+ssgradtube*tubetranene(itype(i)) +!C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- & + ((bordtubetop-positi)/tubebufthick) + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i)=enetube(i)+sstube*tubetranenepep +!C gg_tube_SC(3,i)=gg_tube_SC(3,i) +!C &+ssgradtube*tubetranene(itype(i)) +!C gg_tube(3,i-1)= gg_tube(3,i-1) +!C &+ssgradtube*tubetranene(itype(i)) +!C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i)=enetube(i)+sstube*tubetranenepep +!C print *,"I am in true lipid" + endif + else +!C sstube=0.0d0 +!C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor + +!C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=enetube(i)+sstube* & + (pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6) +!C write(iout,*) "TU13",i,rdiff6,enetube(i) +!C print *,rdiff,rdiff6,pep_aa_tube +!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +!C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6- & + 6.0d0*pep_bb_tube)/rdiff6/rdiff*sstube +!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +!C &rdiff,fac + +!C now direction of gg_tube vector + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + gg_tube(3,i)=gg_tube(3,i) & + +ssgradtube*enetube(i)/sstube/2.0d0 + gg_tube(3,i-1)= gg_tube(3,i-1) & + +ssgradtube*enetube(i)/sstube/2.0d0 + + enddo +!C basically thats all code now we split for side-chains (REMEMBER to sum up at the END) +!C print *,gg_tube(1,0),"TU" + do i=itube_start,itube_end +!C Lets not jump over memory as we use many times iti + iti=itype(i) +!C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) & +!!C in UNRES uncomment the line below as GLY has no side-chain... + .or.(iti.eq.10) & + ) cycle + vectube(1)=c(1,i+nres) + vectube(1)=mod(vectube(1),boxxsize) + if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize + vectube(2)=c(2,i+nres) + vectube(2)=mod(vectube(2),boxysize) + if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) +!C THIS FRAGMENT MAKES TUBE FINITE + 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 +!C respos=mod(c(3,i+nres),boxzsize) +!C print *,positi,bordtubebot,buftubebot,bordtubetop + + if ((positi.gt.bordtubebot) & + .and.(positi.lt.bordtubetop)) then +!C the energy transfer exist + if (positi.lt.buftubebot) then + fracinbuf=1.0d0- & + ((positi-bordtubebot)/tubebufthick) +!C lipbufthick is thickenes of lipid buffore + sstube=sscalelip(fracinbuf) + ssgradtube=-sscagradlip(fracinbuf)/tubebufthick +!C print *,ssgradtube, sstube,tubetranene(itype(i)) + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +!C gg_tube_SC(3,i)=gg_tube_SC(3,i) +!C &+ssgradtube*tubetranene(itype(i)) +!C gg_tube(3,i-1)= gg_tube(3,i-1) +!C &+ssgradtube*tubetranene(itype(i)) +!C print *,"doing sccale for lower part" + elseif (positi.gt.buftubetop) then + fracinbuf=1.0d0- & + ((bordtubetop-positi)/tubebufthick) + + sstube=sscalelip(fracinbuf) + ssgradtube=sscagradlip(fracinbuf)/tubebufthick + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +!C gg_tube_SC(3,i)=gg_tube_SC(3,i) +!C &+ssgradtube*tubetranene(itype(i)) +!C gg_tube(3,i-1)= gg_tube(3,i-1) +!C &+ssgradtube*tubetranene(itype(i)) +!C print *, "doing sscalefor top part",sslip,fracinbuf + else + sstube=1.0d0 + ssgradtube=0.0d0 + enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i)) +!C print *,"I am in true lipid" + endif + else +!C sstube=0.0d0 +!C ssgradtube=0.0d0 + cycle + endif ! if in lipid or buffor +!CEND OF FINITE FRAGMENT +!C as the tube is infinity we do not calculate the Z-vector use of Z +!C as chosen axis + vectube(3)=0.0d0 +!C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +!C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r +!C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +!C and its 6 power + rdiff6=rdiff**6.0d0 +!C for vectorization reasons we will sumup at the end to avoid depenence of previous + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6)& + *sstube+enetube(i+nres) +!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +!C now we calculate gradient + fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff-& + 6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube +!C now direction of gg_tube vector + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + gg_tube_SC(3,i)=gg_tube_SC(3,i) & + +ssgradtube*enetube(i+nres)/sstube + gg_tube(3,i-1)= gg_tube(3,i-1) & + +ssgradtube*enetube(i+nres)/sstube + + enddo + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + enddo +!C print *,"ETUBE", etube + return + end subroutine calctube2 +!===================================================================================================================================== + subroutine calcnano(Etube) + real(kind=8) :: vectube(3),enetube(nres*2), & + enecavtube(nres*2) + real(kind=8) :: Etube,xtemp,xminact,yminact,& + ytemp,xmin,ymin,tub_r,rdiff,rdiff6,fac,denominator,faccav,& + sc_aa_tube,sc_bb_tube,zmin,ztemp,zminact + integer:: i,j,iti + + Etube=0.0d0 + print *,itube_start,itube_end,"poczatek" + do i=itube_start,itube_end + enetube(i)=0.0d0 + enetube(i+nres)=0.0d0 + enddo +!C first we calculate the distance from tube center +!C first sugare-phosphate group for NARES this would be peptide group +!C for UNRES + do i=itube_start,itube_end +!C lets ommit dummy atoms for now + if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle +!C now calculate distance from center of tube and direction vectors + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + + do j=-1,1 + vectube(1)=dmod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=dmod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + vectube(3)=dmod((c(3,i)+c(3,i+1))/2.0d0,boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=dabs(vectube(1)-tubecenter(1)) + yminact=dabs(vectube(2)-tubecenter(2)) + zminact=dabs(vectube(3)-tubecenter(3)) + + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + if (zmin.gt.zminact) then + zmin=zminact + ztemp=vectube(3) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(3)=ztemp + + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + vectube(3)=vectube(3)-tubecenter(3) + +!C print *,"x",(c(1,i)+c(1,i+1))/2.0d0,tubecenter(1) +!C print *,"y",(c(2,i)+c(2,i+1))/2.0d0,tubecenter(2) +!C as the tube is infinity we do not calculate the Z-vector use of Z +!C as chosen axis +!C vectube(3)=0.0d0 +!C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +!C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r + vectube(3)=vectube(3)/tub_r +!C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +!C and its 6 power + rdiff6=rdiff**6.0d0 +!C for vectorization reasons we will sumup at the end to avoid depenence of previous + enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6 +!C write(iout,*) "TU13",i,rdiff6,enetube(i) +!C print *,rdiff,rdiff6,pep_aa_tube +!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +!C now we calculate gradient + fac=(-12.0d0*pep_aa_tube/rdiff6- & + 6.0d0*pep_bb_tube)/rdiff6/rdiff +!C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), +!C &rdiff,fac + if (acavtubpep.eq.0.0d0) then +!C go to 667 + enecavtube(i)=0.0 + faccav=0.0 + else + denominator=(1.0d0+dcavtubpep*rdiff6*rdiff6) + enecavtube(i)= & + (bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff)+ccavtubpep) & + /denominator + enecavtube(i)=0.0 + faccav=((bcavtubpep*1.0d0+acavtubpep/2.0d0/dsqrt(rdiff)) & + *denominator-(bcavtubpep*rdiff+acavtubpep*dsqrt(rdiff) & + +ccavtubpep)*rdiff6**2.0d0/rdiff*dcavtubpep*12.0d0) & + /denominator**2.0d0 +!C faccav=0.0 +!C fac=fac+faccav +!C 667 continue + endif + + do j=1,3 + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac/2.0d0 + gg_tube(j,i)=gg_tube(j,i)+vectube(j)*fac/2.0d0 + enddo + enddo + + do i=itube_start,itube_end + enecavtube(i)=0.0d0 +!C Lets not jump over memory as we use many times iti + iti=itype(i) +!C lets ommit dummy atoms for now + if ((iti.eq.ntyp1) & +!C in UNRES uncomment the line below as GLY has no side-chain... +!C .or.(iti.eq.10) + ) cycle + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + do j=-1,1 + vectube(1)=dmod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=dmod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + vectube(3)=dmod((c(3,i+nres)),boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=dabs(vectube(1)-tubecenter(1)) + yminact=dabs(vectube(2)-tubecenter(2)) + zminact=dabs(vectube(3)-tubecenter(3)) + + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + if (zmin.gt.zminact) then + zmin=zminact + ztemp=vectube(3) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp + vectube(3)=ztemp + +!C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +!C & tubecenter(2) + vectube(1)=vectube(1)-tubecenter(1) + vectube(2)=vectube(2)-tubecenter(2) + vectube(3)=vectube(3)-tubecenter(3) +!C now calculte the distance + tub_r=dsqrt(vectube(1)**2+vectube(2)**2+vectube(3)**2) +!C now normalize vector + vectube(1)=vectube(1)/tub_r + vectube(2)=vectube(2)/tub_r + vectube(3)=vectube(3)/tub_r + +!C calculte rdiffrence between r and r0 + rdiff=tub_r-tubeR0 +!C and its 6 power + rdiff6=rdiff**6.0d0 + sc_aa_tube=sc_aa_tube_par(iti) + sc_bb_tube=sc_bb_tube_par(iti) + enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 +!C enetube(i+nres)=0.0d0 +!C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 +!C now we calculate gradient + fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- & + 6.0d0*sc_bb_tube/rdiff6/rdiff +!C fac=0.0 +!C now direction of gg_tube vector +!C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) + if (acavtub(iti).eq.0.0d0) then +!C go to 667 + enecavtube(i+nres)=0.0d0 + faccav=0.0d0 + else + denominator=(1.0d0+dcavtub(iti)*rdiff6*rdiff6) + enecavtube(i+nres)= & + (bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff)+ccavtub(iti)) & + /denominator +!C enecavtube(i)=0.0 + faccav=((bcavtub(iti)*1.0d0+acavtub(iti)/2.0d0/dsqrt(rdiff)) & + *denominator-(bcavtub(iti)*rdiff+acavtub(iti)*dsqrt(rdiff) & + +ccavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)*12.0d0) & + /denominator**2.0d0 +!C faccav=0.0 + fac=fac+faccav +!C 667 continue + endif +!C print *,"TUT",i,iti,rdiff,rdiff6,acavtub(iti),denominator, +!C & enecavtube(i),faccav +!C print *,"licz=", +!C & (bcavtub(iti)*rdiff+acavtub(iti)*sqrt(rdiff)+ccavtub(iti)) +!C print *,"finene=",enetube(i+nres)+enecavtube(i) + do j=1,3 + gg_tube_SC(j,i)=gg_tube_SC(j,i)+vectube(j)*fac + gg_tube(j,i-1)=gg_tube(j,i-1)+vectube(j)*fac + enddo + enddo + + + + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) & + +enecavtube(i+nres) + enddo +!C print *,"ETUBE", etube + return + end subroutine calcnano + +!=============================================== !-------------------------------------------------------------------------------- !C first for shielding is setting of function of side-chains @@ -18327,6 +18980,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(gshieldc_ll(3,-1:nres)) allocate(gshieldc_loc_ll(3,-1:nres)) allocate(grad_shield(3,-1:nres)) + allocate(gg_tube_sc(3,-1:nres)) + allocate(gg_tube(3,-1:nres)) !(3,maxres) allocate(grad_shield_side(3,50,nres)) allocate(grad_shield_loc(3,50,nres)) diff --git a/source/unres/io.f90 b/source/unres/io.f90 index 7964cdd..40bcb18 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -733,6 +733,7 @@ call reada(weightcard,'WSHIELD',wshield,0.05D0) call reada(weightcard,'LIPSCALE',lipscale,1.0D0) call reada(weightcard,'WLT',wliptran,1.0D0) + call reada(weightcard,'WTUBE',wtube,1.0d0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) call reada(weightcard,'SCAL14',scal14,0.4D0) diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 0dce2e1..91daac1 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -743,7 +743,8 @@ real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,& dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,& - res1,epsijlip + res1,epsijlip,epspeptube,epssctube,sigmapeptube, & + sigmasctube integer :: ichir1,ichir2 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) @@ -2042,6 +2043,10 @@ allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) + allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),& + dcavtub(ntyp1)) + allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),& + tubetranene(ntyp1)) aa_aq(:,:)=0.0D0 bb_aq(:,:)=0.0D0 aa_lip(:,:)=0.0D0 @@ -2049,7 +2054,13 @@ chi(:,:)=0.0D0 sigma(:,:)=0.0D0 r0(:,:)=0.0D0 - + acavtub(:)=0.0d0 + bcavtub(:)=0.0d0 + ccavtub(:)=0.0d0 + dcavtub(:)=0.0d0 + sc_aa_tube_par(:)=0.0d0 + sc_bb_tube_par(:)=0.0d0 + !-------------------------------- do i=2,ntyp @@ -2133,6 +2144,25 @@ endif enddo enddo + write(iout,*) "tube param" + read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, & + ccavtubpep,dcavtubpep,tubetranenepep + sigmapeptube=sigmapeptube**6 + sigeps=dsign(1.0D0,epspeptube) + epspeptube=dabs(epspeptube) + pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2 + pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube + write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep + do i=1,ntyp + read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), & + ccavtub(i),dcavtub(i),tubetranene(i) + sigmasctube=sigmasctube**6 + sigeps=dsign(1.0D0,epssctube) + epssctube=dabs(epssctube) + sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2 + sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube + write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i) + enddo allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2) bad(:,:)=0.0D0 @@ -2943,6 +2973,20 @@ call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) + call readi(controlcard,'TUBEMOD',tubemode,0) + write (iout,*) TUBEmode,"TUBEMODE" + if (TUBEmode.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + buftubetop=bordtubetop-tubebufthick + endif + ! CUTOFFF ON ELECTROSTATICS call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) @@ -3745,6 +3789,9 @@ open (isidep,file=sidename,status='old',action='read') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + ! print *,"Processor",myrank," opened file ISIDEP" ! print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -3772,6 +3819,8 @@ open (isidep,file=sidename,status='old') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& readonly) @@ -3804,6 +3853,9 @@ open (isidep,file=sidename,status='old',readonly) call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',action='read') + #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') -- 1.7.9.5