From: Adam Sieradzan Date: Sun, 6 Aug 2017 13:15:33 +0000 (+0200) Subject: adding ebend_nucl to UCGM+some further reading X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=bbbdc8e18680625d3004f414aec255e9ca7b7353;p=unres4.git adding ebend_nucl to UCGM+some further reading --- diff --git a/PARAM/TiO2_fin.parm b/PARAM/TiO2_fin.parm index a03bc3b..acbdf67 100644 --- a/PARAM/TiO2_fin.parm +++ b/PARAM/TiO2_fin.parm @@ -1,5 +1,5 @@ 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 + 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 @@ -7,7 +7,7 @@ 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.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 diff --git a/PARAM/angle_06292011.parm b/PARAM/angle_06292011.parm new file mode 100644 index 0000000..73126c8 --- /dev/null +++ b/PARAM/angle_06292011.parm @@ -0,0 +1,13 @@ +1 6 0 0 0 0 +1 1 1 1 1 0 +SSS + 8.81610E+00 +-1.01498E+02 + 1.36850E+02 +-1.32224E+02 + 8.75874E+01 +-4.00674E+01 + 1.01373E+01 + + + diff --git a/PARAM/base-GB-Yi-20110420-corr.parm b/PARAM/base-GB-Yi-20110420-corr.parm new file mode 100644 index 0000000..cbb3c99 --- /dev/null +++ b/PARAM/base-GB-Yi-20110420-corr.parm @@ -0,0 +1,23 @@ +2 +3.448E+00 3.234E+00 5.968E-01 6.212E-01 2.786E-01 4.409E-01 -1.058E-01 1.204E-15 -1.065E+00 0.000E+00 2.992E-02 2.992E-02 ! A-A +3.271E+00 3.215E+00 6.358E-01 6.656E-01 3.303E-01 4.624E-01 -2.033E-01 3.560E-02 -1.232E+00 0.000E+00 2.992E-02 3.054E-02 ! A-G +2.851E+00 3.251E+00 5.841E-01 6.552E-01 2.830E-01 4.741E-01 -6.024E-02 -1.435E-01 -1.071E+00 0.000E+00 2.992E-02 5.985E-02 ! A-C +3.448E+00 3.125E+00 6.956E-01 6.481E-01 4.801E-01 5.648E-01 -7.967E-13 1.588E-01 -7.549E-01 0.000E+00 2.992E-02 1.461E-01 ! A-T +3.448E+00 3.125E+00 6.956E-01 6.481E-01 4.801E-01 5.648E-01 -7.967E-13 1.588E-01 -7.549E-01 0.000E+00 2.992E-02 1.461E-01 ! A-U +3.794E+00 3.168E+00 6.424E-01 6.696E-01 3.915E-01 4.515E-01 -7.988E-02 1.809E-16 -1.096E+00 0.000E+00 3.054E-02 3.054E-02 ! G-G +3.332E+00 3.130E+00 6.468E-01 6.698E-01 3.557E-01 4.221E-01 -9.982E-13 -1.599E-01 -7.488E-01 0.000E+00 3.054E-02 5.985E-02 ! G-C +2.906E+00 3.214E+00 6.287E-01 6.616E-01 3.386E-01 4.956E-01 -2.259E-01 7.638E-02 -1.238E+00 0.000E+00 3.054E-02 1.461E-01 ! G-T +2.906E+00 3.214E+00 6.287E-01 6.616E-01 3.386E-01 4.956E-01 -2.259E-01 7.638E-02 -1.238E+00 0.000E+00 3.054E-02 1.461E-01 ! G-U +2.366E+00 3.241E+00 6.939E-01 7.101E-01 3.766E-01 5.576E-01 -3.317E-13 1.149E-13 -1.150E+00 0.000E+00 5.985E-02 5.985E-02 ! C-C +3.068E+00 3.248E+00 5.086E-01 5.928E-01 2.417E-01 4.332E-01 -2.711E-16 -3.275E-01 -6.115E-01 0.000E+00 5.985E-02 1.461E-01 ! C-T +3.068E+00 3.248E+00 5.086E-01 5.928E-01 2.417E-01 4.332E-01 -2.711E-16 -3.275E-01 -6.115E-01 0.000E+00 5.985E-02 1.461E-01 ! C-U +2.864E+00 3.188E+00 6.657E-01 6.873E-01 4.198E-01 4.553E-01 -2.150E-15 5.481E-01 -1.071E+00 0.000E+00 1.461E-01 1.461E-01 ! T-T +2.864E+00 3.188E+00 6.657E-01 6.873E-01 4.198E-01 4.553E-01 -2.150E-15 5.481E-01 -1.071E+00 0.000E+00 1.461E-01 1.461E-01 ! T-U +2.864E+00 3.188E+00 6.657E-01 6.873E-01 4.198E-01 4.553E-01 -2.150E-15 5.481E-01 -1.071E+00 0.000E+00 1.461E-01 1.461E-01 ! U-U + + + +***************************************************************************************************************************************** +*epsi0 sigma0 chi1 chi2 chip1 chip2 Dij Aij Bij Cij d1 d2 * +* Thymine parameters was used for temperarilly instead of Uridine! * +***************************************************************************************************************************************** diff --git a/PARAM/bond_06302011.parm b/PARAM/bond_06302011.parm new file mode 100644 index 0000000..752d958 --- /dev/null +++ b/PARAM/bond_06302011.parm @@ -0,0 +1,6 @@ +5.862 4.825 99.0 33.0 5.0 ! sugar-sugar +4.513 34.909 285.0 95.0 5.4 ! sugar-adenine +4.620 45.924 302.0 100.7 5.0 ! sugar-guanine +3.950 37.644 259.0 86.3 4.0 ! sugar-cytosine +4.187 23.184 259.0 86.3 4.0 | sugar-thymine +3.930 51.542 262.0 87.3 4.0 | sugar-uracil diff --git a/PARAM/rotamers_06292011.parm b/PARAM/rotamers_06292011.parm new file mode 100644 index 0000000..4ecf477 --- /dev/null +++ b/PARAM/rotamers_06292011.parm @@ -0,0 +1,50 @@ +A + -4.09128 + 9.01905 + -1.60446 + -9.63035 + 19.08409 + -7.45394 + 0.26563 +-11.87042 + -7.80643 +G + -4.93113 + 12.26405 + -2.19708 + -9.66999 + 20.64929 + -8.97987 + 0.43400 +-18.01950 + -6.73950 +C + -8.48015 + 32.37047 + -3.84007 +-18.11192 + 36.21928 +-16.10797 + -0.54304 +-24.37115 +-11.89684 +T +-11.89377 + 41.71295 + -1.08909 +-19.52261 + 41.09039 +-19.56908 + 3.45006 +-30.60118 + -0.86021 +U + -3.81240 + 13.34544 + -1.31020 + -8.35478 + 18.00603 + -7.65220 + -1.92437 + -9.48286 + -1.78505 diff --git a/PARAM/torsion_06292011.parm b/PARAM/torsion_06292011.parm new file mode 100644 index 0000000..6f31e33 --- /dev/null +++ b/PARAM/torsion_06292011.parm @@ -0,0 +1,6 @@ +1 +1 1 1 1 1 +3 0 +1 -8.80190E-01 -2.99791E-01 +1 -2.86844E-02 -1.81187E-01 +1 -6.43850E-02 -8.72961E-02 diff --git a/PARAM/torsion_d.parm b/PARAM/torsion_d.parm new file mode 100644 index 0000000..23a8a01 --- /dev/null +++ b/PARAM/torsion_d.parm @@ -0,0 +1,4 @@ +1 +1 1 1 1 1 +1 0 +-1.0 0.0 diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt index fd62a89..dc8c7f1 100644 --- a/source/unres/CMakeLists.txt +++ b/source/unres/CMakeLists.txt @@ -71,8 +71,8 @@ set(UNRES_MD_SRC3 if (Fortran_COMPILER_NAME STREQUAL "ifort") set (CMAKE_Fortran_FLAGS_RELEASE " ") set (CMAKE_Fortran_FLAGS_DEBUG "-O0 -g ") - set(FFLAGS0 "-fpp -c -O3 -ip " ) -# set(FFLAGS0 "-CB -g -ip -fpp" ) +# set(FFLAGS0 "-fpp -c -O3 -ip " ) + set(FFLAGS0 "-CB -g -ip -fpp" ) set(FFLAGS1 "-fpp -c -O " ) set(FFLAGS2 "-fpp -c -g -CA -CB ") set(FFLAGS3 "-fpp -c -g -O0 " ) diff --git a/source/unres/MD.f90 b/source/unres/MD.f90 index 3d3dc39..1cc216e 100644 --- a/source/unres/MD.f90 +++ b/source/unres/MD.f90 @@ -702,7 +702,7 @@ endif ! write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) ! write (iout,*) "i",i," msc",msc(iti)," v",(v(j),j=1,3) - KEt_sc=KEt_sc+msc(iti)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) + KEt_sc=KEt_sc+msc(iti,1)*(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) vtot(i+nres)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) do j=1,3 incr(j)=incr(j)+d_t(j,i) @@ -733,14 +733,14 @@ incr(j)=d_t(j,nres+i) enddo ! write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) - KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ & + KEr_sc=KEr_sc+Isc(iti,1)*(incr(1)*incr(1)+incr(2)*incr(2)+ & incr(3)*incr(3)) endif enddo ! The total kinetic energy 111 continue ! write(iout,*) 'KEr_sc', KEr_sc - KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc) + KE_total=0.5d0*(mp(1)*KEt_p+KEt_sc+0.25d0*Ip(1)*KEr_p+KEr_sc) ! write (iout,*) "KE_total",KE_total return end subroutine kinetic @@ -2881,18 +2881,18 @@ if (i.lt.nct) then do k=1,3 ! write (iout,*) "k",k," ii+k",ii+k," ii+j+k",ii+j+k,"EK1",Ek1 - Ek1=Ek1+0.5d0*mp*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& - 0.5d0*0.25d0*IP*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 + Ek1=Ek1+0.5d0*mp(1)*((d_t_work(ii+j+k)+d_t_work_new(ii+k))/2)**2+& + 0.5d0*0.25d0*IP(1)*(d_t_work(ii+j+k)-d_t_work_new(ii+k))**2 enddo endif if (itype(i,1).ne.10) ii=ii+3 - write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1)) + write (iout,*) "i",i," itype",itype(i,1)," mass",msc(itype(i,1),1) write (iout,*) "ii",ii do k=1,3 ii=ii+1 write (iout,*) "k",k," ii",ii,"EK1",EK1 - if (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1)))*(d_t_work(ii)-d_t_work(ii-3))**2 - Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)))*d_t_work(ii)**2 + if (iabs(itype(i,1)).ne.10) Ek1=Ek1+0.5d0*Isc(iabs(itype(i,1),1))*(d_t_work(ii)-d_t_work(ii-3))**2 + Ek1=Ek1+0.5d0*msc(iabs(itype(i,1)),1)*d_t_work(ii)**2 enddo write (iout,*) "i",i," ii",ii enddo @@ -3678,31 +3678,31 @@ enddo enddo do j=1,3 - cm(j)=mp*cm(j) + cm(j)=mp(1)*cm(j) enddo M_SC=0.0d0 do i=nnt,nct iti=iabs(itype(i,1)) - M_SC=M_SC+msc(iabs(iti)) + M_SC=M_SC+msc(iabs(iti),1) inres=i+nres do j=1,3 - cm(j)=cm(j)+msc(iabs(iti))*c(j,inres) + cm(j)=cm(j)+msc(iabs(iti),1)*c(j,inres) enddo enddo do j=1,3 - cm(j)=cm(j)/(M_SC+(nct-nnt)*mp) + cm(j)=cm(j)/(M_SC+(nct-nnt)*mp(1)) enddo do i=nnt,nct-1 do j=1,3 pr(j)=c(j,i)+0.5d0*dc(j,i)-cm(j) enddo - Im(1,1)=Im(1,1)+mp*(pr(2)*pr(2)+pr(3)*pr(3)) - Im(1,2)=Im(1,2)-mp*pr(1)*pr(2) - Im(1,3)=Im(1,3)-mp*pr(1)*pr(3) - Im(2,3)=Im(2,3)-mp*pr(2)*pr(3) - Im(2,2)=Im(2,2)+mp*(pr(3)*pr(3)+pr(1)*pr(1)) - Im(3,3)=Im(3,3)+mp*(pr(1)*pr(1)+pr(2)*pr(2)) + Im(1,1)=Im(1,1)+mp(1)*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-mp(1)*pr(1)*pr(2) + Im(1,3)=Im(1,3)-mp(1)*pr(1)*pr(3) + Im(2,3)=Im(2,3)-mp(1)*pr(2)*pr(3) + Im(2,2)=Im(2,2)+mp(1)*(pr(3)*pr(3)+pr(1)*pr(1)) + Im(3,3)=Im(3,3)+mp(1)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct @@ -3711,26 +3711,26 @@ do j=1,3 pr(j)=c(j,inres)-cm(j) enddo - Im(1,1)=Im(1,1)+msc(iabs(iti))*(pr(2)*pr(2)+pr(3)*pr(3)) - Im(1,2)=Im(1,2)-msc(iabs(iti))*pr(1)*pr(2) - Im(1,3)=Im(1,3)-msc(iabs(iti))*pr(1)*pr(3) - Im(2,3)=Im(2,3)-msc(iabs(iti))*pr(2)*pr(3) - Im(2,2)=Im(2,2)+msc(iabs(iti))*(pr(3)*pr(3)+pr(1)*pr(1)) - Im(3,3)=Im(3,3)+msc(iabs(iti))*(pr(1)*pr(1)+pr(2)*pr(2)) + Im(1,1)=Im(1,1)+msc(iabs(iti),1)*(pr(2)*pr(2)+pr(3)*pr(3)) + Im(1,2)=Im(1,2)-msc(iabs(iti),1)*pr(1)*pr(2) + Im(1,3)=Im(1,3)-msc(iabs(iti),1)*pr(1)*pr(3) + Im(2,3)=Im(2,3)-msc(iabs(iti),1)*pr(2)*pr(3) + Im(2,2)=Im(2,2)+msc(iabs(iti),1)*(pr(3)*pr(3)+pr(1)*pr(1)) + Im(3,3)=Im(3,3)+msc(iabs(iti),1)*(pr(1)*pr(1)+pr(2)*pr(2)) enddo do i=nnt,nct-1 - Im(1,1)=Im(1,1)+Ip*(1-dc_norm(1,i)*dc_norm(1,i))* & + Im(1,1)=Im(1,1)+Ip(1)*(1-dc_norm(1,i)*dc_norm(1,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(1,2)=Im(1,2)+Ip*(-dc_norm(1,i)*dc_norm(2,i))* & + Im(1,2)=Im(1,2)+Ip(1)*(-dc_norm(1,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(1,3)=Im(1,3)+Ip*(-dc_norm(1,i)*dc_norm(3,i))* & + Im(1,3)=Im(1,3)+Ip(1)*(-dc_norm(1,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(2,3)=Im(2,3)+Ip*(-dc_norm(2,i)*dc_norm(3,i))* & + Im(2,3)=Im(2,3)+Ip(1)*(-dc_norm(2,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(2,2)=Im(2,2)+Ip*(1-dc_norm(2,i)*dc_norm(2,i))* & + Im(2,2)=Im(2,2)+Ip(1)*(1-dc_norm(2,i)*dc_norm(2,i))* & vbld(i+1)*vbld(i+1)*0.25d0 - Im(3,3)=Im(3,3)+Ip*(1-dc_norm(3,i)*dc_norm(3,i))* & + Im(3,3)=Im(3,3)+Ip(1)*(1-dc_norm(3,i)*dc_norm(3,i))* & vbld(i+1)*vbld(i+1)*0.25d0 enddo @@ -3739,17 +3739,17 @@ if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then iti=iabs(itype(i,1)) inres=i+nres - Im(1,1)=Im(1,1)+Isc(iti)*(1-dc_norm(1,inres)* & + Im(1,1)=Im(1,1)+Isc(iti,1)*(1-dc_norm(1,inres)* & dc_norm(1,inres))*vbld(inres)*vbld(inres) - Im(1,2)=Im(1,2)-Isc(iti)*(dc_norm(1,inres)* & + Im(1,2)=Im(1,2)-Isc(iti,1)*(dc_norm(1,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) - Im(1,3)=Im(1,3)-Isc(iti)*(dc_norm(1,inres)* & + Im(1,3)=Im(1,3)-Isc(iti,1)*(dc_norm(1,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) - Im(2,3)=Im(2,3)-Isc(iti)*(dc_norm(2,inres)* & + Im(2,3)=Im(2,3)-Isc(iti,1)*(dc_norm(2,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) - Im(2,2)=Im(2,2)+Isc(iti)*(1-dc_norm(2,inres)* & + Im(2,2)=Im(2,2)+Isc(iti,1)*(1-dc_norm(2,inres)* & dc_norm(2,inres))*vbld(inres)*vbld(inres) - Im(3,3)=Im(3,3)+Isc(iti)*(1-dc_norm(3,inres)* & + Im(3,3)=Im(3,3)+Isc(iti,1)*(1-dc_norm(3,inres)* & dc_norm(3,inres))*vbld(inres)*vbld(inres) endif enddo @@ -3871,7 +3871,7 @@ enddo call vecpr(pr(1),v(1),vp) do j=1,3 - L(j)=L(j)+mp*vp(j) + L(j)=L(j)+mp(1)*vp(j) enddo do j=1,3 pr(j)=0.5d0*dc(j,i) @@ -3879,7 +3879,7 @@ enddo call vecpr(pr(1),pp(1),vp) do j=1,3 - L(j)=L(j)+Ip*vp(j) + L(j)=L(j)+Ip(1)*vp(j) enddo enddo do j=1,3 @@ -3904,7 +3904,7 @@ ! write (iout,*) "i",i," iti",iti," pr",(pr(j),j=1,3), ! & " v",(v(j),j=1,3)," vp",(vp(j),j=1,3) do j=1,3 - L(j)=L(j)+msc(iabs(iti))*vp(j) + L(j)=L(j)+msc(iabs(iti),1)*vp(j) enddo ! write (iout,*) "L",(l(j),j=1,3) if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then @@ -3913,7 +3913,7 @@ enddo call vecpr(dc(1,inres),d_t(1,inres),vp) do j=1,3 - L(j)=L(j)+Isc(iti)*vp(j) + L(j)=L(j)+Isc(iti,1)*vp(j) enddo endif do j=1,3 @@ -3947,12 +3947,12 @@ summas=0.0d0 do i=nnt,nct if (i.lt.nct) then - summas=summas+mp + summas=summas+mp(1) do j=1,3 - vcm(j)=vcm(j)+mp*(vv(j)+0.5d0*d_t(j,i)) + vcm(j)=vcm(j)+mp(1)*(vv(j)+0.5d0*d_t(j,i)) enddo endif - amas=msc(iabs(itype(i,1))) + amas=msc(iabs(itype(i,1)),1) summas=summas+amas if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then do j=1,3 @@ -5361,12 +5361,12 @@ enddo ! Load peptide group radii do i=nnt,nct-1 - radius(i)=pstok + radius(i)=pstok(1) enddo ! Load side chain radii do i=nnt,nct iti=itype(i,1) - radius(i+nres)=restok(iti) + radius(i+nres)=restok(iti,1) enddo ! do i=1,2*nres ! write (iout,*) "i",i," radius",radius(i) diff --git a/source/unres/REMD.f90 b/source/unres/REMD.f90 index 6a558a4..8a049a7 100644 --- a/source/unres/REMD.f90 +++ b/source/unres/REMD.f90 @@ -502,8 +502,8 @@ do i=nnt,nct-1 ind=ind+1 ind1=ind1+1 - coeff=0.25d0*IP - massvec(ind1)=mp + coeff=0.25d0*IP(1) + massvec(ind1)=mp(1) Gmat(ind,ind)=coeff A(ind1,ind)=0.5d0 enddo @@ -520,17 +520,17 @@ m1=nct-nnt+1 ind=0 ind1=0 - msc(ntyp1)=1.0d0 + msc(ntyp1,1)=1.0d0 do i=nnt,nct ind=ind+1 ii = ind+m iti=itype(i,1) - massvec(ii)=msc(iabs(iti)) + massvec(ii)=msc(iabs(iti),1) if (iti.ne.10 .and. iti.ne.ntyp1) then ind1=ind1+1 ii1= ind1+m1 A(ii,ii1)=1.0d0 - Gmat(ii1,ii1)=ISC(iabs(iti)) + Gmat(ii1,ii1)=ISC(iabs(iti),1) endif enddo ! Off-diagonal elements of the dX part of A diff --git a/source/unres/compare.F90 b/source/unres/compare.F90 index 24d2aa6..34f8a2a 100644 --- a/source/unres/compare.F90 +++ b/source/unres/compare.F90 @@ -471,7 +471,7 @@ integer :: ncont,i,j,i1,j1,nbeta,nstrand,ii1,jj1,ij,nhelix,& iii1,jjj1 integer,dimension(2,12*nres) :: icont !(2,maxcont) (maxcont=12*maxres) - integer,dimension(nres,4) :: isec !(maxres,4) + integer,dimension(nres,0:4) :: isec !(maxres,4) integer,dimension(nres) :: nsec !(maxres) logical :: lprint,not_done !,freeres real(kind=8) :: p1,p2 diff --git a/source/unres/control.F90 b/source/unres/control.F90 index 38ac803..a45ba29 100644 --- a/source/unres/control.F90 +++ b/source/unres/control.F90 @@ -229,6 +229,16 @@ iweight=31 izsc=32 #endif + ibond_nucl=126 + ithep_nucl=127 + irotam_nucl=128 + itorp_nucl= 129 + itordp_nucl= 130 +! ielep_nucl= 131 + isidep_nucl=132 + + + iliptranpar=60 itube=61 #if defined(WHAM_RUN) || defined(CLUSTER) @@ -477,7 +487,7 @@ logical :: scheck,lprint,flag !el local variables - integer :: ind_scint=0,ind_scint_old,ii,jj,i,j,iint + integer :: ind_scint=0,ind_scint_old,ii,jj,i,j,iint,itmp #ifdef MPI integer :: my_sc_int(0:nfgtasks-1),my_ele_int(0:nfgtasks-1) @@ -490,7 +500,7 @@ ind_scpint_old,nsumgrad,nlen,ngrad_start,ngrad_end,& ierror,k,ierr,iaux,ncheck_to,ncheck_from,ind_typ,& ichunk,int_index_old - integer,dimension(5) :: nct_molec + integer,dimension(5) :: nct_molec,nnt_molec !el allocate(itask_cont_from(0:nfgtasks-1)) !(0:max_fg_procs-1) !el allocate(itask_cont_to(0:nfgtasks-1)) !(0:max_fg_procs-1) @@ -502,10 +512,33 @@ itask_cont_to(i)=fg_rank enddo lprint=energy_dec + itmp=0 + do i=1,5 + if (nres_molec(i).eq.0) cycle + itmp=itmp+nres_molec(i) + if (itype(itmp,i).eq.ntyp1_molec(i)) then + nct_molec(i)=itmp-1 + else + nct_molec(i)=itmp + endif + enddo +! nct_molec(1)=nres_molec(1)-1 + itmp=0 + do i=2,5 + itmp=itmp+nres_molec(i-1) + if (itype(itmp+1,i).eq.ntyp1_molec(i)) then + nnt_molec(i)=itmp+2 + else + nnt_molec(i)=itmp+1 + endif + enddo + print *,"nres_molec",nres_molec(:) + print *,"nnt_molec",nnt_molec(:) + print *,"nct_molec",nct_molec(:) ! lprint=.true. if (lprint) & write (iout,*)'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct - n_sc_int_tot=(nct-nnt+1)*(nct-nnt)/2-nss + n_sc_int_tot=(nct_molec(1)-nnt+1)*(nct_molec(1)-nnt)/2-nss call int_bounds(n_sc_int_tot,my_sc_inds,my_sc_inde) !write (iout,*) 'INIT_INT_TABLE nres=',nres,' nnt=',nnt,' nct=',nct if (lprint) & @@ -533,6 +566,7 @@ allocate(itask_cont_to_all(0:nfgtasks-1,0:nfgtasks-1)) !el---------- ! lprint=.false. + print *,"NCT",nct_molec(1),nct do i=1,nres !el !maxres nint_gr(i)=0 nscp_gr(i)=0 @@ -549,7 +583,9 @@ ind_scint_old=0 !d write (iout,*) 'ns=',ns,' nss=',nss,' ihpb,jhpb', !d & (ihpb(i),jhpb(i),i=1,nss) - do i=nnt,nct-1 +! print *,nnt,nct_molec(1) + do i=nnt,nct_molec(1)-1 +! print*, "inloop",i scheck=.false. if (dyn_ss) goto 10 do ii=1,nss @@ -560,27 +596,28 @@ endif enddo 10 continue +! print *,'i=',i,' scheck=',scheck,' jj=',jj !d write (iout,*) 'i=',i,' scheck=',scheck,' jj=',jj if (scheck) then if (jj.eq.i+1) then #ifdef MPI ! write (iout,*) 'jj=i+1' call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,& - iatsc_s,iatsc_e,i+2,nct,nint_gr(i),istart(i,1),iend(i,1),*12) + iatsc_s,iatsc_e,i+2,nct_molec(1),nint_gr(i),istart(i,1),iend(i,1),*12) #else nint_gr(i)=1 istart(i,1)=i+2 iend(i,1)=nct #endif - else if (jj.eq.nct) then + else if (jj.eq.nct_molec(1)) then #ifdef MPI ! write (iout,*) 'jj=nct' call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,& - iatsc_s,iatsc_e,i+1,nct-1,nint_gr(i),istart(i,1),iend(i,1),*12) + iatsc_s,iatsc_e,i+1,nct_molec(1)-1,nint_gr(i),istart(i,1),iend(i,1),*12) #else nint_gr(i)=1 istart(i,1)=i+1 - iend(i,1)=nct-1 + iend(i,1)=nct_molecule(1)-1 #endif else #ifdef MPI @@ -589,6 +626,7 @@ ii=nint_gr(i)+1 call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,& iatsc_s,iatsc_e,jj+1,nct_molec(1),nint_gr(i),istart(i,ii),iend(i,ii),*12) + #else nint_gr(i)=2 istart(i,1)=i+1 @@ -599,8 +637,12 @@ endif else #ifdef MPI +! print *,"i for EVDW",iatsc_s,iatsc_e,istart(i,1),iend(i,1),& +! i+1,nct_molec(1),nint_gr(i),ind_scint,my_sc_inds,my_sc_inde,i call int_partition(ind_scint,my_sc_inds,my_sc_inde,i,& - iatsc_s,iatsc_e,i+1,nct_molec(1),nint_gr(i),istart(i,1),iend(i,1),*12) + iatsc_s,iatsc_e,i+1,nct_molec(1),nint_gr(i), & + istart(i,1),iend(i,1),*12) +! print *,"i for EVDW",iatsc_s,iatsc_e,istart(i,1),iend(i,1) #else nint_gr(i)=1 istart(i,1)=i+1 @@ -613,6 +655,8 @@ #endif enddo 12 continue +! print *,"i for EVDW",iatsc_s,iatsc_e,istart(i,1),iend(i,1) + #ifndef MPI iatsc_s=nnt iatsc_e=nct-1 @@ -622,6 +666,7 @@ if (lprint) write (*,*) 'Processor',fg_rank,' CG Group',kolor,& ' absolute rank',myrank,' iatsc_s=',iatsc_s,' iatsc_e=',iatsc_e #endif +! lprint=.true. if (lprint) then write (iout,'(a)') 'Interaction array:' do i=iatsc_s,iatsc_e @@ -629,6 +674,7 @@ i,(istart(i,iint),iend(i,iint),iint=1,nint_gr(i)) enddo endif +! lprint=.false. ispp=4 !?? wham ispp=2 #ifdef MPI ! Now partition the electrostatic-interaction array @@ -648,11 +694,11 @@ iatel_e=0 ind_eleint=0 ind_eleint_old=0 - if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then - nct_molec(1)=nres_molec(1)-1 - else - nct_molec(1)=nres_molec(1) - endif +! if (itype(nres_molec(1),1).eq.ntyp1_molec(1)) then +! nct_molec(1)=nres_molec(1)-1 +! else +! nct_molec(1)=nres_molec(1) +! endif ! print *,"nct",nct,nct_molec(1),itype(nres_molec(1),1),ntyp_molec(1) do i=nnt,nct_molec(1)-3 ijunk=0 @@ -776,6 +822,9 @@ call int_bounds(nres_molec(1)-2,ithet_start,ithet_end) ithet_start=ithet_start+2 ithet_end=ithet_end+2 + call int_bounds(nres_molec(2)-2,ithet_nucl_start,ithet_nucl_end) + ithet_nucl_start=ithet_nucl_start+2+nres_molec(1) + ithet_nucl_end=ithet_nucl_end+2+nres_molec(1) call int_bounds(nct_molec(1)-nnt-2,iturn3_start,iturn3_end) iturn3_start=iturn3_start+nnt iphi_start=iturn3_start+2 @@ -800,10 +849,21 @@ call int_bounds(nres_molec(1)-2,ibond_start,ibond_end) ibond_start=ibond_start+1 ibond_end=ibond_end+1 - print *,ibond_start,ibond_end +! print *,ibond_start,ibond_end call int_bounds(nct_molec(1)-nnt,ibondp_start,ibondp_end) ibondp_start=ibondp_start+nnt ibondp_end=ibondp_end+nnt + call int_bounds(nres_molec(2)-2,ibond_nucl_start,ibond_nucl_end) + ibond_nucl_start=ibond_nucl_start+nnt_molec(2)-1 + ibond_nucl_end=ibond_nucl_end+nnt_molec(2)-1 + print *,"NUCLibond",ibond_nucl_start,ibond_nucl_end + print *, "before devision",nnt_molec(2),nct_molec(2)-nnt_molec(2) + call int_bounds(nct_molec(2)-nnt_molec(2),ibondp_nucl_start,ibondp_nucl_end) + ibondp_nucl_start=ibondp_nucl_start+nnt_molec(2) + ibondp_nucl_end=ibondp_nucl_end+nnt_molec(2) + print *,"NUCLibond2",ibondp_nucl_start,ibondp_nucl_end + + call int_bounds1(nres_molec(1)-1,ivec_start,ivec_end) ! print *,"Processor",myrank,fg_rank,fg_rank1, ! & " ivec_start",ivec_start," ivec_end",ivec_end @@ -1324,6 +1384,8 @@ loc_end=nres_molec(1)-1 ithet_start=3 ithet_end=nres_molec(1) + ithet_nucl_start=3+nres_molec(1) + ithet_nucl_end=nres_molec(1)+nres_molec(2) iturn3_start=nnt iturn3_end=nct_molec(1)-3 iturn4_start=nnt @@ -1342,8 +1404,12 @@ itau_end=nres_molec(1) ibond_start=2 ibond_end=nres_molec(1)-1 + ibond_nucl_start=2+nres_molec(1) + ibond_nucl_end=nres_molec(2)-1 ibondp_start=nnt ibondp_end=nct_molec(1)-1 + ibondp_nucl_start=nnt_molec(2) + ibondp_nucl_end=nct_molec(2) ivec_start=1 ivec_end=nres_molec(1)-1 iset_start=3 @@ -1747,7 +1813,7 @@ else if (molecule.eq.2) then do i=1,ntyp1_molec(molecule) print *,nam(1:1),restyp(i,molecule)(1:1) - if (nam(1:1).eq.restyp(i,molecule)(1:1)) then + if (nam(2:2).eq.restyp(i,molecule)(1:1)) then rescode=i return endif @@ -1760,8 +1826,8 @@ stop else if (molecule.eq.5) then do i=1,ntyp1_molec(molecule) - print *,i,restyp(i,molecule) - if (ucase(nam).eq.restyp(i,molecule)) then + print *,i,restyp(i,molecule)(1:2) + if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then rescode=i return endif diff --git a/source/unres/data/MD_data.f90 b/source/unres/data/MD_data.f90 index 2bfe6e5..6087b9a 100644 --- a/source/unres/data/MD_data.f90 +++ b/source/unres/data/MD_data.f90 @@ -64,15 +64,16 @@ integer :: dimen,dimen1,dimen3 integer :: lang,count_reset_moment,count_reset_vel ! common /inertia/ - real(kind=8) :: IP,mp - real(kind=8),dimension(:),allocatable :: ISC,msc !(ntyp+1) + real(kind=8),dimension(5) :: IP,mp + real(kind=8),dimension(:,:),allocatable :: ISC,msc !(ntyp+1) ! common /langevin/ - real(kind=8) :: rwat,etawat,stdfp,pstok,gamp!,Rb + real(kind=8) :: rwat,etawat,stdfp,gamp!,Rb + real(kind=8), dimension(5) :: pstok real(kind=8) :: cPoise=2.9361d0, Rb=0.001986d0 real(kind=8),dimension(:),allocatable :: gamsc !(ntyp1) real(kind=8),dimension(:),allocatable :: stdfsc !(ntyp) - real(kind=8),dimension(:),allocatable :: restok !(ntyp+1) + real(kind=8),dimension(:,:),allocatable :: restok !(ntyp+1) logical :: surfarea integer :: reset_fricmat ! common /mdpmpi/ diff --git a/source/unres/data/energy_data.f90 b/source/unres/data/energy_data.f90 index a0e0f3a..c41fc11 100644 --- a/source/unres/data/energy_data.f90 +++ b/source/unres/data/energy_data.f90 @@ -64,7 +64,8 @@ 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,wtube + wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, & + wbond_nucl,wang_nucl #ifdef CLUSTER real(kind=8) :: scalscp #endif @@ -100,7 +101,7 @@ integer,dimension(:,:),allocatable :: istart,iend !(maxres,maxint_gr) integer,dimension(:),allocatable :: nint_gr,itel,& ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr !(maxres) - integer,dimension(:),allocatable :: istype + integer,dimension(:),allocatable :: istype,molnum integer,dimension(:,:),allocatable :: itype ! now itype has more molecule types integer,dimension(:,:),allocatable :: iscpstart,iscpend !(maxres,maxint_gr) integer :: iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,& @@ -118,7 +119,13 @@ ! common /stretch/ real(kind=8) :: vbldp0,akp,distchainmax,vbldpDUM real(kind=8),dimension(:,:),allocatable :: vbldsc0,aksc,abond0 !(maxbondterm,ntyp) + real(kind=8) :: vbldp0_nucl,akp_nucl + real(kind=8),dimension(:,:),allocatable :: vbldsc0_nucl,& + aksc_nucl,abond0_nucl !(maxbondterm,ntyp) + integer,dimension(:),allocatable :: nbondterm !(ntyp) + integer,dimension(:),allocatable :: nbondterm_nucl !(ntyp) + !----------------------------------------------------------------------------- ! common.local ! Parameters of ab initio-derived potential of virtual-bond-angle bending @@ -133,6 +140,20 @@ ccthet,ddthet,eethet !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) real(kind=8),dimension(:,:,:,:,:,:,:),allocatable :: ffthet,ggthet + +!-----------nucleic acid parameters-------------------------- + integer :: nthetyp_nucl,ntheterm_nucl,ntheterm2_nucl,& + ntheterm3_nucl,nsingle_nucl,& + ndouble_nucl,nntheterm_nucl + integer,dimension(:),allocatable :: ithetyp_nucl !(-ntyp1:ntyp1) + real(kind=8),dimension(:,:,:),allocatable :: aa0thet_nucl +!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + real(kind=8),dimension(:,:,:,:),allocatable :: aathet_nucl + real(kind=8),dimension(:,:,:,:,:),allocatable :: bbthet_nucl,& + ccthet_nucl,ddthet_nucl,eethet_nucl +!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + real(kind=8),dimension(:,:,:,:,:,:),allocatable :: ffthet_nucl,ggthet_nucl + !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) ! Parameters of the virtual-bond-angle probability distribution ! common /thetas/ @@ -158,6 +179,8 @@ iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start,& iint_end,iphi1_start,iphi1_end,itau_start,itau_end,& ilip_start,ilip_end,itube_start,itube_end + integer :: ibond_nucl_start,ibond_nucl_end, & + ibondp_nucl_start,ibondp_nucl_end,ithet_nucl_start,ithet_nucl_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 ebda347..8b79b35 100644 --- a/source/unres/data/io_units.f90 +++ b/source/unres/data/io_units.f90 @@ -14,7 +14,9 @@ integer :: inp,iout,igeom,intin,ipdb,imol2,ipdbin,ithep,irotam,& itorp,itordp,ifourier,ielep,isidep,iscpp,icbase,istat,ientin,& ientout,izs1,isecpred,ibond,irest2,iifrag,icart,irest1,isccor,& - ithep_pdb,irotam_pdb,iliptranpar,itube + ithep_pdb,irotam_pdb,iliptranpar,itube, & + ibond_nucl,ithep_nucl,irotam_nucl,itorp_nucl,itordp_nucl, & + isidep_nucl #ifdef WHAM_RUN ! el wham iounits integer :: isidep1,ihist,iweight,izsc,idistr @@ -43,7 +45,9 @@ ! common /parfiles/ character(len=256) :: bondname,thetname,rotname,torname,tordname,& fouriername,elename,sidename,scpname,sccorname,patname,& - thetname_pdb,rotname_pdb,liptranname,tubename + thetname_pdb,rotname_pdb,liptranname,tubename, & + bondname_nucl,thetname_nucl,rotname_nucl,torname_nucl, & + tordname_nucl,sidename_nucl !----------------------------------------------------------------------- ! INP - main input file ! IOUT - list file diff --git a/source/unres/data/names.f90 b/source/unres/data/names.f90 index e669cc1..97eb2d0 100644 --- a/source/unres/data/names.f90 +++ b/source/unres/data/names.f90 @@ -59,7 +59,7 @@ !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Number of energy components - integer,parameter :: n_ene=25 + integer,parameter :: n_ene=38 integer :: n_ene2=2*n_ene !----------------------------------------------------------------------------- ! common.names @@ -83,17 +83,24 @@ "ECORR6 ","EELLO ","ETURN3 ","ETURN4 ","ETURN6 ",& "EBE bend ","ESC SCloc ","ETORS ","ETORSD ","EHPB ","EVDWPP ",& "EVDW2_14 ","ESTR ","ESCCOR ","EDIHC ","EVDW_T ",& - "ELT "," "," ","ETUBE " /) + "ELT "," "," ","ETUBE ",& + "EVDWPP ","EESPP ","EVDWPSB ","EESPSB ","EVDWSB ",& + "EESSB ","ESTR ","EBE ","ESBLOC ","ETORS ",& + "ETORSD ","ECORR ","ECORR3 " /) 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 ",& - "WLT "," "," ","WTUBE " /) + "WLT "," "," ","WTUBE " ,& + "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",& + "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",& + "WTORD ","WCORR ","WCORR3 "/) integer :: nprint_ene = 21 integer,dimension(n_ene) :: print_order = & - (/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25/) + (/1,2,3,18,11,12,13,14,4,5,6,7,8,9,10,19,16,15,17,20,21,22,23,24,25,& + 26,27,28,29,30,31,32,33,34,35,36,37,38/) !#endif diff --git a/source/unres/energy.f90 b/source/unres/energy.f90 index 3975993..2fc2773 100644 --- a/source/unres/energy.f90 +++ b/source/unres/energy.f90 @@ -125,6 +125,8 @@ gshieldc_t3,gshieldc_loc_t3,gshieldx_t4,gshieldc_t4, & gshieldc_loc_t4,gshieldx_ll,gshieldc_ll,gshieldc_loc_ll,& grad_shield,gg_tube,gg_tube_sc,gradafm !(3,maxres) +!-----------------------------NUCLEIC GRADIENT + real(kind=8),dimension(:,:),allocatable ::gradb_nucl,gradbx_nucl ! 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) @@ -226,7 +228,10 @@ real(kind=8) :: eello_turn3,eello_turn4,estr,ebe,eliptran,etube, & Eafmforce,ethetacnstr real(kind=8) :: ecorr,ecorr5,ecorr6,eturn6 - +! now energies for nulceic alone parameters + real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& + ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& + ecorr3_nucl #ifdef MPI real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw ! shielding effect varibles for MPI @@ -357,6 +362,7 @@ if (shield_mode.eq.2) then call set_shield_fac2 endif + print *,"AFTER EGB",ipot,evdw !mc !mc Sep-06: egb takes care of dynamic ss bonds too !mc @@ -539,7 +545,10 @@ else etube=0.0d0 endif - +!-------------------------------------------------------- + call ebond_nucl(estr_nucl) + call ebend_nucl(ebe_nucl) + print *,"after ebend", ebe_nucl #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -585,6 +594,21 @@ energia(23)=Eafmforce energia(24)=ethetacnstr energia(25)=etube +!--------------------------------------------------------------- + energia(26)=evdwpp + energia(27)=eespp + energia(28)=evdwpsb + energia(29)=eelpsb + energia(30)=evdwsb + energia(31)=eelsb + energia(32)=estr_nucl + energia(33)=ebe_nucl + energia(34)=esbloc + energia(35)=etors_nucl + energia(36)=etors_d_nucl + energia(37)=ecorr_nucl + energia(38)=ecorr3_nucl +!---------------------------------------------------------------------- ! 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" @@ -627,6 +651,10 @@ real(kind=8) :: eel_loc,eello_turn3,eello_turn4,eturn6,ebe,escloc real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,esccor,etot, & eliptran,etube, Eafmforce,ethetacnstr + real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& + ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& + ecorr3_nucl + integer :: i #ifdef MPI integer :: ierr @@ -690,6 +718,9 @@ Eafmforce=energia(23) ethetacnstr=energia(24) etube=energia(25) + estr_nucl=energia(32) + ebe_nucl=energia(33) + #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & @@ -697,7 +728,8 @@ +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+wtube*etube& - +Eafmforce+ethetacnstr + +Eafmforce+ethetacnstr & + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & @@ -705,8 +737,8 @@ +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+wtube*etube& - +Eafmforce+ethetacnstr - + +Eafmforce+ethetacnstr & + +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl #endif energia(0)=etot ! detecting NaNQ @@ -831,6 +863,9 @@ real(kind=8) :: eello_turn6,eello_turn3,eello_turn4,ebe,escloc real(kind=8) :: etors,etors_d,ehpb,edihcnstr,estr,Uconst,esccor,eliptran,& etube,ethetacnstr,Eafmforce + real(kind=8) :: evdwpp,eespp,evdwpsb,eelpsb,evdwsb,eelsb,estr_nucl,& + ebe_nucl,esbloc,etors_nucl,etors_d_nucl,ecorr_nucl,& + ecorr3_nucl etot=energia(0) evdw=energia(1) @@ -864,6 +899,9 @@ Eafmforce=energia(23) ethetacnstr=energia(24) etube=energia(25) + estr_nucl=energia(32) + ebe_nucl=energia(33) + #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,& estr,wbond,ebe,wang,& @@ -872,7 +910,9 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,& edihcnstr,ethetacnstr,ebr*nss,& - Uconst,eliptran,wliptran,Eafmforce,etube,wtube,etot + Uconst,eliptran,wliptran,Eafmforce,etube,wtube, & ! till now protein + estr_nucl,wbond_nucl,ebe_nucl,wang_nucl, & + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -900,6 +940,8 @@ 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/& 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & + 'ESTR_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ & + 'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ & 'ETOT= ',1pE16.6,' (total)') #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,& @@ -909,7 +951,9 @@ ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,edihcnstr,& ethetacnstr,ebr*nss,Uconst,eliptran,wliptran,Eafmforc, & - etube,wtube,etot + etube,wtube, & + estr_nucl,wbond_nucl, ebe_nucl,wang_nucl,& + etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ & @@ -936,6 +980,8 @@ 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETUBE=',1pE16.6, ' WEIGHT=',1pD16.6,' (cylindrical energy)'/ & + 'ESTR_nucl= ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching for nucleic)'/ & + 'EBE_nucl=',1pE16.6,' WEIGHT=',1pD16.6,' (bending for nucleic)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1610,6 +1656,7 @@ !C print *,i,j,c(1,i),c(1,j),c(2,i),c(2,j),c(3,i),c(3,j) ! if (energy_dec) write (iout,*) & ! 'evdw',i,j,evdwij +! print *,"ZALAMKA", evdw ! Calculate gradient components. e1=e1*eps1*eps2rt**2*eps3rt**2 @@ -1642,6 +1689,7 @@ enddo ! j enddo ! iint enddo ! i +! print *,"ZALAMKA", evdw ! write (iout,*) "Number of loop steps in EGB:",ind !ccc energy_dec=.false. return @@ -5145,7 +5193,7 @@ ! endif enddo estr=0.5d0*AKP*estr+estr1 - print *,"estr_bb",estr,AKP +! print *,"estr_bb",estr,AKP ! ! 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included ! @@ -5160,7 +5208,7 @@ "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,& AKSC(1,iti),AKSC(1,iti)*diff*diff estr=estr+0.5d0*AKSC(1,iti)*diff*diff - print *,"estr_sc",estr +! print *,"estr_sc",estr do j=1,3 gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) enddo @@ -5189,11 +5237,11 @@ usumsqder=usumsqder+ud(j)*uprod2 enddo estr=estr+uprod/usum - print *,"estr_sc",estr,i +! print *,"estr_sc",estr,i if (energy_dec) write (iout,*) & "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,& - AKSC(1,iti),AKSC(1,iti)*diff*diff + AKSC(1,iti),uprod/usum do j=1,3 gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) enddo @@ -10290,10 +10338,8 @@ +wturn3*gshieldc_t3(j,i)& +wturn4*gshieldc_t4(j,i)& +wel_loc*gshieldc_ll(j,i)& - +wtube*gg_tube(j,i) - - - + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) enddo enddo #else @@ -10315,9 +10361,8 @@ +wcorr*gshieldc_ec(j,i) & +wturn4*gshieldc_t4(j,i) & +wel_loc*gshieldc_ll(j,i)& - +wtube*gg_tube(j,i) - - + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) enddo enddo @@ -10474,7 +10519,9 @@ +wturn4*gshieldc_loc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) & +wel_loc*gshieldc_loc_ll(j,i) & - +wtube*gg_tube(j,i) + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) + #else @@ -10508,7 +10555,9 @@ +wturn4*gshieldc_loc_t4(j,i) & +wel_loc*gshieldc_ll(j,i) & +wel_loc*gshieldc_loc_ll(j,i) & - +wtube*gg_tube(j,i) + +wtube*gg_tube(j,i) & + +wbond_nucl*gradb_nucl(j,i) + @@ -10524,7 +10573,9 @@ +wturn3*gshieldx_t3(j,i) & +wturn4*gshieldx_t4(j,i) & +wel_loc*gshieldx_ll(j,i)& - +wtube*gg_tube_sc(j,i) + +wtube*gg_tube_sc(j,i) & + +wbond_nucl*gradbx_nucl(j,i) + enddo @@ -10734,7 +10785,7 @@ ! include 'COMMON.CALC' ! include 'COMMON.IOUNITS' real(kind=8), dimension(3) :: dcosom1,dcosom2 - +! print *,"wchodze" eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1 eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2 eom12=evdwij*eps1_om12+eps2der*eps2rt_om12 & @@ -16106,6 +16157,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' gg_tube(j,i)=0.0d0 gg_tube_sc(j,i)=0.0d0 gradafm(j,i)=0.0d0 + gradb_nucl(j,i)=0.0d0 + gradbx_nucl(j,i)=0.0d0 do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo @@ -18740,7 +18793,7 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' 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 + integer:: i,j,iti,r Etube=0.0d0 ! print *,itube_start,itube_end,"poczatek" @@ -18948,6 +19001,23 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' Etube=Etube+enetube(i)+enetube(i+nres)+enecavtube(i) & +enecavtube(i+nres) enddo +! do i=1,20 +! print *,"begin", i,"a" +! do r=1,10000 +! rdiff=r/100.0d0 +! rdiff6=rdiff**6.0d0 +! sc_aa_tube=sc_aa_tube_par(i) +! sc_bb_tube=sc_bb_tube_par(i) +! enetube(i)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 +! denominator=(1.0d0+dcavtub(i)*rdiff6*rdiff6) +! enecavtube(i)= & +! (bcavtub(i)*rdiff+acavtub(i)*dsqrt(rdiff)+ccavtub(i)) & +! /denominator + +! print '(5(f10.3,1x))',rdiff,enetube(i),enecavtube(i),enecavtube(i)+enetube(i) +! enddo +! print *,"end",i,"a" +! enddo !C print *,"ETUBE", etube return end subroutine calcnano @@ -19435,6 +19505,8 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' allocate(gg_tube_sc(3,-1:nres)) allocate(gg_tube(3,-1:nres)) allocate(gradafm(3,-1:nres)) + allocate(gradb_nucl(3,-1:nres)) + allocate(gradbx_nucl(3,-1:nres)) !(3,maxres) allocate(grad_shield_side(3,50,nres)) allocate(grad_shield_loc(3,50,nres)) @@ -19612,6 +19684,283 @@ write(iout,*) 'Calling CHECK_ECARTIN else.' return end subroutine alloc_ener_arrays +!----------------------------------------------------------------- + subroutine ebond_nucl(estr_nucl) +!c +!c Evaluate the energy of stretching of the CA-CA and CA-SC virtual bonds +!c + + real(kind=8),dimension(3) :: u,ud + real(kind=8) :: usum,uprod,uprod1,uprod2,usumsqder + real(kind=8) :: estr_nucl,diff + integer :: iti,i,j,k,nbi + estr_nucl=0.0d0 +!C print *,"I enter ebond" + if (energy_dec) & + write (iout,*) "ibondp_start,ibondp_end",& + ibondp_nucl_start,ibondp_nucl_end + do i=ibondp_nucl_start,ibondp_nucl_end + if (itype(i-1,2).eq.ntyp1_molec(2) .or. & + itype(i,2).eq.ntyp1_molec(2)) cycle +! estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax) +! do j=1,3 +! gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) +! & *dc(j,i-1)/vbld(i) +! enddo +! if (energy_dec) write(iout,*) +! & "estr1",i,vbld(i),distchainmax, +! & gnmr1(vbld(i),-1.0d0,distchainmax) + + diff = vbld(i)-vbldp0_nucl + if(energy_dec)write(iout,*) "estr_nucl_bb" , i,vbld(i),& + vbldp0_nucl,diff,AKP_nucl*diff*diff + estr_nucl=estr_nucl+diff*diff + print *,estr_nucl + do j=1,3 + gradb_nucl(j,i-1)=AKP_nucl*diff*dc(j,i-1)/vbld(i) + enddo +!c write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3) + enddo + estr_nucl=0.5d0*AKP_nucl*estr_nucl + print *,"partial sum", estr_nucl,AKP_nucl + + if (energy_dec) & + write (iout,*) "ibondp_start,ibondp_end",& + ibond_nucl_start,ibond_nucl_end + + do i=ibond_nucl_start,ibond_nucl_end +!C print *, "I am stuck",i + iti=itype(i,2) + if (iti.eq.ntyp1_molec(2)) cycle + nbi=nbondterm_nucl(iti) +!C print *,iti,nbi + if (nbi.eq.1) then + diff=vbld(i+nres)-vbldsc0_nucl(1,iti) + + if (energy_dec) & + write (iout,*) "estr_nucl_sc", i,iti,vbld(i+nres),vbldsc0_nucl(1,iti),diff, & + AKSC_nucl(1,iti),AKSC_nucl(1,iti)*diff*diff + estr_nucl=estr_nucl+0.5d0*AKSC_nucl(1,iti)*diff*diff + print *,estr_nucl + do j=1,3 + gradbx_nucl(j,i)=AKSC_nucl(1,iti)*diff*dc(j,i+nres)/vbld(i+nres) + enddo + else + do j=1,nbi + diff=vbld(i+nres)-vbldsc0_nucl(j,iti) + ud(j)=aksc_nucl(j,iti)*diff + u(j)=abond0_nucl(j,iti)+0.5d0*ud(j)*diff + enddo + uprod=u(1) + do j=2,nbi + uprod=uprod*u(j) + enddo + usum=0.0d0 + usumsqder=0.0d0 + do j=1,nbi + uprod1=1.0d0 + uprod2=1.0d0 + do k=1,nbi + if (k.ne.j) then + uprod1=uprod1*u(k) + uprod2=uprod2*u(k)*u(k) + endif + enddo + usum=usum+uprod1 + usumsqder=usumsqder+ud(j)*uprod2 + enddo + estr_nucl=estr_nucl+uprod/usum + do j=1,3 + gradbx(j,i)=usumsqder/(usum*usum)*dc(j,i+nres)/vbld(i+nres) + enddo + endif + enddo +!C print *,"I am about to leave ebond" + return + end subroutine ebond_nucl + +!----------------------------------------------------------------------------- + subroutine ebend_nucl(etheta_nucl) + real(kind=8),dimension(nntheterm_nucl+1) :: coskt,sinkt !mmaxtheterm + real(kind=8),dimension(nsingle_nucl+1) :: cosph1,sinph1,cosph2,sinph2 !maxsingle + real(kind=8),dimension(ndouble_nucl+1,ndouble_nucl+1) :: cosph1ph2,sinph1ph2 !maxdouble,maxdouble + logical :: lprn=.true., lprn1=.false. +!el local variables + integer :: i,k,iblock,ityp1,ityp2,ityp3,l,m + real(kind=8) :: dethetai,dephii,dephii1,theti2,phii,phii1,ethetai + real(kind=8) :: aux,etheta_nucl,ccl,ssl,scl,csl,ethetacnstr +! local variables for constrains + real(kind=8) :: difi,thetiii + integer itheta + etheta_nucl=0.0D0 + print *,"ithet_start",ithet_nucl_start," ithet_end",ithet_nucl_end,nres + do i=ithet_nucl_start,ithet_nucl_end + if ((itype(i-1,2).eq.ntyp1_molec(2)).or.& + (itype(i-2,2).eq.ntyp1_molec(2)).or. & + (itype(i,2).eq.ntyp1_molec(2))) cycle + dethetai=0.0d0 + dephii=0.0d0 + dephii1=0.0d0 + theti2=0.5d0*theta(i) + ityp2=ithetyp_nucl(itype(i-1,2)) + do k=1,nntheterm_nucl + coskt(k)=dcos(k*theti2) + sinkt(k)=dsin(k*theti2) + enddo + if (i.gt.3 .and. itype(i-2,2).ne.ntyp1_molec(2)) then +#ifdef OSF + phii=phi(i) + if (phii.ne.phii) phii=150.0 +#else + phii=phi(i) +#endif + ityp1=ithetyp_nucl(itype(i-2,2)) + do k=1,nsingle_nucl + cosph1(k)=dcos(k*phii) + sinph1(k)=dsin(k*phii) + enddo + else + phii=0.0d0 + ityp1=nthetyp_nucl+1 + do k=1,nsingle_nucl + cosph1(k)=0.0d0 + sinph1(k)=0.0d0 + enddo + endif + + if (i.lt.nres .and. itype(i,2).ne.ntyp1_molec(2)) then +#ifdef OSF + phii1=phi(i+1) + if (phii1.ne.phii1) phii1=150.0 + phii1=pinorm(phii1) +#else + phii1=phi(i+1) +#endif + ityp3=ithetyp_nucl(itype(i,2)) + do k=1,nsingle_nucl + cosph2(k)=dcos(k*phii1) + sinph2(k)=dsin(k*phii1) + enddo + else + phii1=0.0d0 + ityp3=nthetyp_nucl+1 + do k=1,nsingle_nucl + cosph2(k)=0.0d0 + sinph2(k)=0.0d0 + enddo + endif + ethetai=aa0thet_nucl(ityp1,ityp2,ityp3) + do k=1,ndouble_nucl + do l=1,k-1 + ccl=cosph1(l)*cosph2(k-l) + ssl=sinph1(l)*sinph2(k-l) + scl=sinph1(l)*cosph2(k-l) + csl=cosph1(l)*sinph2(k-l) + cosph1ph2(l,k)=ccl-ssl + cosph1ph2(k,l)=ccl+ssl + sinph1ph2(l,k)=scl+csl + sinph1ph2(k,l)=scl-csl + enddo + enddo + if (lprn) then + write (iout,*) "i",i," ityp1",ityp1," ityp2",ityp2,& + " ityp3",ityp3," theti2",theti2," phii",phii," phii1",phii1 + write (iout,*) "coskt and sinkt",nntheterm_nucl + do k=1,nntheterm_nucl + write (iout,*) k,coskt(k),sinkt(k) + enddo + endif + do k=1,ntheterm_nucl + ethetai=ethetai+aathet_nucl(k,ityp1,ityp2,ityp3)*sinkt(k) + dethetai=dethetai+0.5d0*k*aathet_nucl(k,ityp1,ityp2,ityp3)& + *coskt(k) + if (lprn)& + write (iout,*) "k",k," aathet",aathet_nucl(k,ityp1,ityp2,ityp3),& + " ethetai",ethetai + enddo + if (lprn) then + write (iout,*) "cosph and sinph" + do k=1,nsingle_nucl + write (iout,*) k,cosph1(k),sinph1(k),cosph2(k),sinph2(k) + enddo + write (iout,*) "cosph1ph2 and sinph2ph2" + do k=2,ndouble_nucl + do l=1,k-1 + write (iout,*) l,k,cosph1ph2(l,k),cosph1ph2(k,l),& + sinph1ph2(l,k),sinph1ph2(k,l) + enddo + enddo + write(iout,*) "ethetai",ethetai + endif + do m=1,ntheterm2_nucl + do k=1,nsingle_nucl + aux=bbthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)& + +ccthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)& + +ddthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)& + +eethet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k) + ethetai=ethetai+sinkt(m)*aux + dethetai=dethetai+0.5d0*m*aux*coskt(m) + dephii=dephii+k*sinkt(m)*(& + ccthet_nucl(k,m,ityp1,ityp2,ityp3)*cosph1(k)-& + bbthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph1(k)) + dephii1=dephii1+k*sinkt(m)*(& + eethet_nucl(k,m,ityp1,ityp2,ityp3)*cosph2(k)-& + ddthet_nucl(k,m,ityp1,ityp2,ityp3)*sinph2(k)) + if (lprn) & + write (iout,*) "m",m," k",k," bbthet",& + bbthet_nucl(k,m,ityp1,ityp2,ityp3)," ccthet",& + ccthet_nucl(k,m,ityp1,ityp2,ityp3)," ddthet",& + ddthet_nucl(k,m,ityp1,ityp2,ityp3)," eethet",& + eethet_nucl(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai + enddo + enddo + if (lprn) & + write(iout,*) "ethetai",ethetai + do m=1,ntheterm3_nucl + do k=2,ndouble_nucl + do l=1,k-1 + aux=ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+& + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l) + ethetai=ethetai+sinkt(m)*aux + dethetai=dethetai+0.5d0*m*coskt(m)*aux + dephii=dephii+l*sinkt(m)*(& + -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-& + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + dephii1=dephii1+(k-l)*sinkt(m)*( & + -ffthet_nucl(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+& + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)) + if (lprn) then + write (iout,*) "m",m," k",k," l",l," ffthet", & + ffthet_nucl(l,k,m,ityp1,ityp2,ityp3), & + ffthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ggthet",& + ggthet_nucl(l,k,m,ityp1,ityp2,ityp3),& + ggthet_nucl(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai + write (iout,*) cosph1ph2(l,k)*sinkt(m), & + cosph1ph2(k,l)*sinkt(m),& + sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m) + endif + enddo + enddo + enddo +10 continue + if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') & + i,theta(i)*rad2deg,phii*rad2deg, & + phii1*rad2deg,ethetai + etheta_nucl=etheta_nucl+ethetai + print *,i,"partial sum",etheta_nucl + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang_nucl*dephii + if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang_nucl*dephii1 + gloc(nphi+i-2,icg)=wang_nucl*dethetai + enddo + return + end subroutine ebend_nucl + !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module energy diff --git a/source/unres/geometry.f90 b/source/unres/geometry.f90 index 8542f52..6895472 100644 --- a/source/unres/geometry.f90 +++ b/source/unres/geometry.f90 @@ -2851,6 +2851,7 @@ endif endif do i=1,nres-1 + if (molnum(i).ne.1) cycle !in wham do i=1,nres iti=itype(i,1) if ((dist(i,i+1).lt.2.0D0 .or. dist(i,i+1).gt.5.0D0).and.& @@ -2903,8 +2904,9 @@ di=dist(i,nres+i) !#ifndef WHAM_RUN ! 10/03/12 Adam: Correction for zero SC-SC bond length + if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1 .and. di.eq.0.0d0) & - di=dsc(itype(i,1)) + di=dsc(itype(i,molnum(i))) vbld(i+nres)=di if (itype(i,1).ne.10) then vbld_inv(i+nres)=1.0d0/di @@ -2916,12 +2918,21 @@ alph(i)=alpha(nres+i,i,nres2+2) omeg(i)=beta(nres+i,i,nres2+2,i+1) endif + if (iti.ne.0) then if(me.eq.king.or..not.out1file)then if (lprn) & write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),& rad2deg*theta(i),rad2deg*phi(i),dsc(iti),vbld(nres+i),& rad2deg*alph(i),rad2deg*omeg(i) endif + else + if(me.eq.king.or..not.out1file)then + if (lprn) & + write (iout,'(a3,i4,7f10.3)') restyp(iti,1),i,vbld(i),& + rad2deg*theta(i),rad2deg*phi(i),dsc(iti+1),vbld(nres+i),& + rad2deg*alph(i),rad2deg*omeg(i) + endif + endif enddo else if (lprn) then do i=2,nres @@ -3043,7 +3054,8 @@ do j=1,3 sccmj=0.0D0 do i=1,nscat - sccmj=sccmj+sccor(j,i) + sccmj=sccmj+sccor(j,i) + print *,"insccent", ires,sccor(j,i) enddo dc(j,ires)=sccmj/nscat enddo diff --git a/source/unres/io.f90 b/source/unres/io.f90 index ae6e636..9388f15 100644 --- a/source/unres/io.f90 +++ b/source/unres/io.f90 @@ -784,6 +784,7 @@ call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) if (index(weightcard,'SOFT').gt.0) ipot=6 + call reada(weightcard,'WBOND_NUCL',wbond_nucl,1.0D0) ! 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) if (wcorr4.gt.0.0d0) wcorr=wcorr4 diff --git a/source/unres/io_config.f90 b/source/unres/io_config.f90 index 5153cc9..25a5e24 100644 --- a/source/unres/io_config.f90 +++ b/source/unres/io_config.f90 @@ -776,12 +776,15 @@ ! allocate(dsc(ntyp1)) !(ntyp1) allocate(dsc_inv(ntyp1)) !(ntyp1) + allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp) + allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) + allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) allocate(nbondterm(ntyp)) !(ntyp) allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp) - allocate(msc(ntyp+1)) !(ntyp+1) - allocate(isc(ntyp+1)) !(ntyp+1) - allocate(restok(ntyp+1)) !(ntyp+1) + allocate(msc(ntyp+1,5)) !(ntyp+1) + allocate(isc(ntyp+1,5)) !(ntyp+1) + allocate(restok(ntyp+1,5)) !(ntyp+1) allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(long_r_sidechain(ntyp)) allocate(short_r_sidechain(ntyp)) @@ -801,10 +804,10 @@ endif enddo #else - read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok - do i=1,ntyp + read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1) + do i=1,ntyp_molec(1) read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),& - j=1,nbondterm(i)),msc(i),isc(i),restok(i) + j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1) dsc(i) = vbldsc0(1,i) if (i.eq.10) then dsc_inv(i)=0.0D0 @@ -813,14 +816,38 @@ endif enddo #endif + read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2) + do i=1,ntyp_molec(2) + nbondterm_nucl(i)=1 + read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2) +! dsc(i) = vbldsc0_nucl(1,i) +! if (i.eq.10) then +! dsc_inv(i)=0.0D0 +! else +! dsc_inv(i)=1.0D0/dsc(i) +! endif + enddo +! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2) +! do i=1,ntyp_molec(2) +! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),& +! aksc_nucl(j,i),abond0_nucl(j,i),& +! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2) +! dsc(i) = vbldsc0(1,i) +! if (i.eq.10) then +! dsc_inv(i)=0.0D0 +! else +! dsc_inv(i)=1.0D0/dsc(i) +! endif +! enddo + if (lprint) then write(iout,'(/a/)')"Dynamic constants of the interaction sites:" write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',& 'inertia','Pstok' - write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok + write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1) do i=1,ntyp write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),& - vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i) + vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1) do j=2,nbondterm(i) write (iout,'(13x,3f10.5)') & vbldsc0(j,i),aksc(j,i),abond0(j,i) @@ -1207,6 +1234,65 @@ close (ithep_pdb) #endif close(ithep) +!--------------- Reading theta parameters for nucleic acid------- + read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,& + ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl + nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl) + allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) + allocate(aa0thet_nucl(nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(maxtheterm,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1)) +!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) + allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) + allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& + nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) + +!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& +! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) + + read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2)) + + aa0thet_nucl(:,:,:)=0.0d0 + aathet_nucl(:,:,:,:)=0.0d0 + bbthet_nucl(:,:,:,:,:)=0.0d0 + ccthet_nucl(:,:,:,:,:)=0.0d0 + ddthet_nucl(:,:,:,:,:)=0.0d0 + eethet_nucl(:,:,:,:,:)=0.0d0 + ffthet_nucl(:,:,:,:,:,:)=0.0d0 + ggthet_nucl(:,:,:,:,:,:)=0.0d0 + + do i=1,nthetyp_nucl + do j=1,nthetyp_nucl + do k=1,nthetyp_nucl + read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3 + read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k) + read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl) + read (ithep_nucl,*,end=111,err=111) & + (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & + (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl) + read (ithep_nucl,*,end=111,err=111) & + (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), & + ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), & + llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl) + enddo + enddo + enddo !------------------------------------------- allocate(nlob(ntyp1)) !(ntyp1) @@ -2337,7 +2423,7 @@ character(len=80) :: card real(kind=8),dimension(3,20) :: sccor integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, - integer :: isugar + integer :: isugar,molecprev character*1 :: sugar real(kind=8) :: cou real(kind=8),dimension(3) :: ccc @@ -2476,7 +2562,7 @@ ! nres_molec(molecule)=nres_molec(molecule)+1 else molecule=2 - itype(ires,molecule)=rescode(ires,res(3:4),0,molecule) + itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) ! nres_molec(molecule)=nres_molec(molecule)+1 endif endif @@ -2546,7 +2632,8 @@ atom.ne.'OXT' .and. atom(:2).ne.'3H' & .and. atom.ne.'P '.and. & atom(1:1).ne.'H' .and. & - atom.ne.'OP1' .and. atom.ne.'OP2 ') then + atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'& + ) then ! write (iout,*) "sidechain ",atom ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) if ((molecule.ne.2).or.(atom(3:3).ne."'")) then @@ -2555,7 +2642,23 @@ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) endif endif - endif + else if ((ions).and.(card(1:6).eq.'HETATM')) then + + read (card(12:16),*) atom + print *,"HETATOM", atom + read (card(18:20),'(a3)') res + if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.& + (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') & + .or.(atom(1:2).eq.'K ')) & + then + ires=ires+1 + if (molecule.ne.5) molecprev=molecule + molecule=5 + nres_molec(molecule)=nres_molec(molecule)+1 + itype(ires,molecule)=rescode(ires,res(2:4),0,molecule) + read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) + endif + endif !atom enddo 10 write (iout,'(a,i5)') ' Number of residues found: ',ires if (ires.eq.0) return @@ -2565,7 +2668,7 @@ if (ires_old.ne.ires) nres_molec(molecule)=nres_molec(molecule)-2 ! print *,'I have', nres_molec(:) - do k=1,5 + do k=1,4 ! ions are without dummy if (nres_molec(k).eq.0) cycle do i=2,nres-1 ! write (iout,*) i,itype(i,1) @@ -2655,8 +2758,9 @@ nsup=nres nstart_sup=1 ! print *,"molecule",molecule - if (itype(nres,1).ne.10) then + if ((itype(nres,1).ne.10)) then nres=nres+1 + if (molecule.eq.5) molecule=molecprev itype(nres,molecule)=ntyp1_molec(molecule) nres_molec(molecule)=nres_molec(molecule)+1 if (unres_pdb) then @@ -2800,6 +2904,7 @@ ! NOW LETS ROCK! SORTING allocate(c_temporary(3,2*nres)) allocate(itype_temporary(nres,5)) + allocate(molnum(nres)) itype_temporary(:,:)=0 seqalingbegin=1 do k=1,5 @@ -2811,6 +2916,7 @@ enddo itype_temporary(seqalingbegin,k)=itype(i,k) + molnum(seqalingbegin)=k seqalingbegin=seqalingbegin+1 endif enddo @@ -3544,11 +3650,11 @@ if(me.eq.king.or..not.out1file) & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",& eta - gamp=scal_fric*(pstok+rwat)*eta + gamp=scal_fric*(pstok(i)+rwat)*eta stdfp=dsqrt(2*Rb*t_bath/d_time) allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1) do i=1,ntyp - gamsc(i)=scal_fric*(restok(i)+rwat)*eta + gamsc(i)=scal_fric*(restok(i,1)+rwat)*eta stdfsc(i)=dsqrt(2*Rb*t_bath/d_time) enddo if(me.eq.king.or..not.out1file)then @@ -3557,7 +3663,7 @@ " stochastic forces of fully exposed sites" write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp) do i=1,ntyp - write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i),& + write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& gamsc(i),stdfsc(i)*dsqrt(gamsc(i)) enddo endif @@ -3919,6 +4025,8 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',readonly,shared) + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared) call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',readonly,shared) call getenv_loc('ROTPAR',rotname) @@ -3933,6 +4041,19 @@ open (ielep,file=elename,status='old',readonly,shared) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly,shared) + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared) + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared) + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',readonly,shared) + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared) + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared) + + #elif (defined CRAY) || (defined AIX) open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',& action='read') @@ -3943,6 +4064,9 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',action='read') + ! print *,"Processor",myrank," opened file IBOND" call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') @@ -3967,6 +4091,18 @@ ! print *,"Processor",myrank," opened file IELEP" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',action='read') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',action='read') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',action='read') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) @@ -3981,6 +4117,9 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old') + call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old') call getenv_loc('ROTPAR',rotname) @@ -3997,6 +4136,17 @@ open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + + open (ithep_nucl,file=thetname_nucl,status='old') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') call getenv_loc('TUBEPAR',tubename) @@ -4009,6 +4159,8 @@ ! Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) open (ibond,file=bondname,status='old',action='read') + call getenv_loc('BONDPAR_NUCL',bondname_nucl) + open (ibond_nucl,file=bondname_nucl,status='old',action='read') call getenv_loc('THETPAR',thetname) open (ithep,file=thetname,status='old',action='read') call getenv_loc('ROTPAR',rotname) @@ -4031,6 +4183,18 @@ open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) + + call getenv_loc('THETPAR_NUCL',thetname_nucl) + open (ithep_nucl,file=thetname_nucl,status='old',action='read') + call getenv_loc('ROTPAR_NUCL',rotname_nucl) + open (irotam_nucl,file=rotname_nucl,status='old',action='read') + call getenv_loc('TORPAR_NUCL',torname_nucl) + open (itorp_nucl,file=torname_nucl,status='old',action='read') + call getenv_loc('TORDPAR_NUCL',tordname_nucl) + open (itordp_nucl,file=tordname_nucl,status='old',action='read') + call getenv_loc('SIDEPAR_NUCL',sidename_nucl) + open (isidep_nucl,file=sidename_nucl,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) diff --git a/source/unres/md_calc.f90 b/source/unres/md_calc.f90 index 50f23d7..8b8e04d 100644 --- a/source/unres/md_calc.f90 +++ b/source/unres/md_calc.f90 @@ -2213,7 +2213,7 @@ Q = 0 ! ********** ESTABLISH AND PROCESS NEXT SUBMATRIX ********** 100 P = Q + 1 - IP = P + 1 + IP(1) = P + 1 ! DO 120 Q = P, N IF (Q .EQ. N) GO TO 140 @@ -2236,7 +2236,7 @@ GO TO 600 160 NORM = ABS(D(P)) ! - DO 180 I = IP, Q + DO 180 I = IP(1), Q 180 NORM = NORM + ABS(D(I)) + ABS(E(I)) ! ********** EPS2 IS THE CRITERION FOR GROUPING, ! EPS3 REPLACES ZERO PIVOTS AND EQUAL @@ -2327,7 +2327,7 @@ 460 RV6(I) = RV6(I) * XU ! ********** ELIMINATION OPERATIONS ON NEXT VECTOR ! ITERATE ********** - 480 DO 520 I = IP, Q + 480 DO 520 I = IP(1), Q U = RV6(I) ! ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE ! WAS PERFORMED EARLIER IN THE