X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMD.f90;h=1cc216e681bd0ad3149ea7adfa9d61585753d03e;hb=bbbdc8e18680625d3004f414aec255e9ca7b7353;hp=3d3dc39ec738a7635562d24521a674cb9963ae72;hpb=7c0faf2ccb6f94ed9a77aa527d3885ba054d3fb2;p=unres4.git 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)