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
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
--- /dev/null
+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
+
+
+
--- /dev/null
+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! *
+*****************************************************************************************************************************************
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+1
+1 1 1 1 1
+1 0
+-1.0 0.0
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 " )
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)
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
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
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
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
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
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)
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
! 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
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
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
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)
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
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
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
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)
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)
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)
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) &
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
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
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
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
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
#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
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
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
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
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
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
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
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
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
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
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/
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
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,&
! 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
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/
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,&
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
! 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
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
! Number of energy components
- integer,parameter :: n_ene=25
+ integer,parameter :: n_ene=38
integer :: n_ene2=2*n_ene
!-----------------------------------------------------------------------------
! common.names
"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
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)
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
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
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
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"
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
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 &
+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 &
+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
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)
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,&
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)'/ &
'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,&
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)'/ &
'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
!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
enddo ! j
enddo ! iint
enddo ! i
+! print *,"ZALAMKA", evdw
! write (iout,*) "Number of loop steps in EGB:",ind
!ccc energy_dec=.false.
return
! 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
!
"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
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
+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
+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
+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
+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)
+
+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
! 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 &
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
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"
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
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))
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
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.&
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
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
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
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
!
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))
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
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)
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)
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
! 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
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
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
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)
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
! 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
enddo
itype_temporary(seqalingbegin,k)=itype(i,k)
+ molnum(seqalingbegin)=k
seqalingbegin=seqalingbegin+1
endif
enddo
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
" 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
! 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)
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')
! 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')
! 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)
! 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)
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)
! 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)
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)
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
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
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
projects / unres4.git / commitdiff