From: Adam Sieradzan Date: Mon, 23 Jan 2017 10:57:16 +0000 (+0100) Subject: introduction nanosphere part 1 X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=d2ee52eab91e6055d9cd36a42447313d73b00ea1 introduction nanosphere part 1 --- diff --git a/PARAM/TiO2.parm b/PARAM/TiO2.parm new file mode 100644 index 0000000..e236309 --- /dev/null +++ b/PARAM/TiO2.parm @@ -0,0 +1 @@ + 1 0.286 4.508 1.491 -0.238 -2.215 3.922E-14 ! Cys 2 4.173 4.136 -143.035 28.494 180.705 1.086E-10 ! Met 1 10.793 4.152 -250.420 48.235 328.091 7.237E-11 ! Phe 2 9.872 3.799 -308.062 62.468 382.714 3.262E-10 ! Ile 2 12.025 4.056 -448.084 90.596 559.212 2.482E-10 ! Leu 1 3.929 3.894 -241.329 50.844 288.017 6.165E-10 ! Val 2 0.285 5.797 738.577 -155.331 -884.672 3.662E-09 ! Trp 1 29.509 4.046 -648.094 126.324 840.706 1.297E-10 ! Tyr 1 0.058 4.710 0.612 -0.108 -0.790 -3.640E-14 ! Ala 0 0.000 0.000 0.000 0.000 0.000 0.000E+00 ! Gly 1 163.545 1.976 -5996.349 1676.896 5422.473 5.557E-13 ! Thr 1 158.959 1.810 -8243.740 2491.600 6900.467 1.632E-11 ! Ser 3 1.582 3.601 673.244 -173.206 -646.238 4.709E-08 ! Gln 1 7.354 3.256 -485.919 112.979 530.772 9.510E-16 ! Asn 3 25.379 2.776 -729.290 168.012 800.982 2.274E-16 ! Glu 1 25.969 2.505 -711.633 167.698 761.605 8.831E-16 ! Asp 1 6.267 3.952 -176.750 35.322 222.723 1.641E-10 ! His 3 0.755 4.649 6.176 -0.945 -9.823 1.249E-13 ! Arg 3 0.001 8.312 2.122 0.132 -6.640 3.808E-11 ! Lys 1 0.000 0.000 0.000 0.000 0.000 0.000E+00 ! Pro \ No newline at end of file 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 38efafa..a675f1f 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -396,6 +396,8 @@ C print *,"just before call" call calctube(Etube) elseif (TUBElog.eq.2) then call calctube2(Etube) + elseif (TUBElog.eq.3) then + call calcnano(Etube) else Etube=0.0d0 endif @@ -672,25 +674,25 @@ c enddo enddo enddo - j=1 - i=0 - print *,"KUPA2",gradbufc(j,i),wsc*gvdwc(j,i), - & wscp*gvdwc_scp(j,i),gvdwc_scpp(j,i), - & welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i), - & wel_loc*gel_loc_long(j,i), - & wcorr*gradcorr_long(j,i), - & wcorr5*gradcorr5_long(j,i), - & wcorr6*gradcorr6_long(j,i), - & wturn6*gcorr6_turn_long(j,i), - & wstrain*ghpbc(j,i) - & ,wliptran*gliptranc(j,i) - & ,gradafm(j,i) - & ,welec*gshieldc(j,i) - & ,wcorr*gshieldc_ec(j,i) - & ,wturn3*gshieldc_t3(j,i) - & ,wturn4*gshieldc_t4(j,i) - & ,wel_loc*gshieldc_ll(j,i) - & ,wtube*gg_tube(j,i) +C j=1 +C i=0 +C print *,"KUPA2",gradbufc(j,i),wsc*gvdwc(j,i), +C & wscp*gvdwc_scp(j,i),gvdwc_scpp(j,i), +C & welec*gelc_long(j,i),wvdwpp*gvdwpp(j,i), +C & wel_loc*gel_loc_long(j,i), +C & wcorr*gradcorr_long(j,i), +C & wcorr5*gradcorr5_long(j,i), +C & wcorr6*gradcorr6_long(j,i), +C & wturn6*gcorr6_turn_long(j,i), +C & wstrain*ghpbc(j,i) +C & ,wliptran*gliptranc(j,i) +C & ,gradafm(j,i) +C & ,welec*gshieldc(j,i) +C & ,wcorr*gshieldc_ec(j,i) +C & ,wturn3*gshieldc_t3(j,i) +C & ,wturn4*gshieldc_t4(j,i) +C & ,wel_loc*gshieldc_ll(j,i) +C & ,wtube*gg_tube(j,i) #else do i=0,nct do j=1,3 @@ -12484,3 +12486,218 @@ C 4) add reading the center of tube C 5) add COMMONs C 6) add to zerograd + +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 calcnano(Etube) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + double precision tub_r,vectube(3),enetube(maxres*2) + 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 + xmin=boxxsize + ymin=boxysize + zmin=boxzsize + + 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 + vectube(3)=mod((c(3,i)+c(3,i+1))/2.0d0,boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + zminact=abs(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 + +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 + + do i=itube_start,itube_end + enecavtube(i)=0.0 +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)=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 + vectube(3)=mod((c(3,i+nres)),boxzsize) + vectube(3)=vectube(3)+boxzsize*j + + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + zminact=abs(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 +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 +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 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) go to 667 + denominator=(1.0+dcavtub(iti)*rdiff6*rdiff6) + enecavtube(i)= + & acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff)+cavtub(iti)) + & /denominator + faccav=(acavtub(iti)*(1.0+bcavtub(iti)/2.0/sqrt(rdiff)) + & *denominator-acavtub(iti)*(rdiff+bcavtub(iti)*sqrt(rdiff) + & +cavtub(iti))*rdiff6**2.0d0/rdiff*dcavtub(iti)) + & /denominator**2.0d0 + fac=fac+faccav + 667 continue + 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 +C Now cavity term E=a(x+bsqrt(x)+c)/(1+dx^12) +C do i=itube_start,itube_end +C enecav(i)=0.0 +C iti=itype(i) +C if (acavtub(iti).eq.0.0) cycle + + + + do i=itube_start,itube_end + Etube=Etube+enetube(i)+enetube(i+nres) + enddo +C print *,"ETUBE", etube + return + end +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 + diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 577c1ee..ce4c59f 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -270,6 +270,7 @@ C endif if (TUBElog.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) @@ -2007,6 +2008,7 @@ c---------------------------------------------------------------------------- iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 + read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end