+#ifdef DEBUG
+ write (iout,*) "Random_vel, fivediag"
+#endif
+ d_t=0.0d0
+ Ek2=0.0d0
+ EK=0.0d0
+ Ek3=0.0d0
+ do ichain=1,nchain
+ ind=0
+ ghalf=0.0d0
+ n=dimen_chain(ichain)
+ innt=iposd_chain(ichain)
+ inct=innt+n-1
+#ifdef DEBUG
+ write (iout,*) "Chain",ichain," n",n," start",innt
+ do i=innt,inct
+ if (i.lt.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i),
+ & DU2orig(i)
+ else if (i.eq.inct-1) then
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i),DU1orig(i)
+ else
+ write (iout,'(2i3,3f10.5)') i,i-innt+1,DMorig(i)
+ endif
+ enddo
+#endif
+ ghalf(ind+1)=dmorig(innt)
+ ghalf(ind+2)=du1orig(innt)
+ ghalf(ind+3)=dmorig(innt+1)
+ ind=ind+3
+ do i=3,n
+ ind=ind+i-3
+c write (iout,*) "i",i," ind",ind," indu2",innt+i-2,
+c & " indu1",innt+i-1," indm",innt+i
+ ghalf(ind+1)=du2orig(innt-1+i-2)
+ ghalf(ind+2)=du1orig(innt-1+i-1)
+ ghalf(ind+3)=dmorig(innt-1+i)
+c write (iout,'(3(a,i2,1x))') "DU2",innt-1+i-2,
+c & "DU1",innt-1+i-1,"DM ",innt-1+i
+ ind=ind+3
+ enddo
+#ifdef DEBUG
+ ind=0
+ do i=1,n
+ do j=1,i
+ ind=ind+1
+ inertia(i,j)=ghalf(ind)
+ inertia(j,i)=ghalf(ind)
+ enddo
+ enddo
+#endif
+#ifdef DEBUG
+ write (iout,*) "Chain ",ichain," ind",ind," dim",n*(n+1)/2
+ write (iout,*) "Five-diagonal inertia matrix, lower triangle"
+ call matoutr(n,ghalf)
+#endif
+ call gldiag(maxres2_chain,n,n,Ghalf,work,Geigen,Gvec,ierr,iwork)
+ if (large) then
+ write (iout,'(//a,i3)')
+ & "Eigenvectors and eigenvalues of the G matrix chain",ichain
+ call eigout(n,n,maxres2_chain,maxres2_chain,Gvec,Geigen)
+ endif
+#ifdef DIAGCHECK
+c check diagonalization
+ do i=1,n
+ do j=1,n
+ aux=0.0d0
+ do k=1,n
+ do l=1,n
+ aux=aux+gvec(k,i)*gvec(l,j)*inertia(k,l)
+ enddo
+ enddo
+ if (i.eq.j) then
+ write (iout,*) i,j,aux,geigen(i)
+ else
+ write (iout,*) i,j,aux
+ endif
+ enddo
+ enddo
+#endif
+ xv=0.0d0
+ ii=0
+ do i=1,n
+ do k=1,3
+ ii=ii+1
+ sigv=dsqrt((Rb*t_bath)/geigen(i))
+ lowb=-5*sigv
+ highb=5*sigv
+ d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+ EK=EK+0.5d0*geigen(i)*d_t_work_new(ii)**2
+c write (iout,*) "i",i," ii",ii," geigen",geigen(i),
+c & " d_t_work_new",d_t_work_new(ii)
+ enddo
+ enddo
+ do k=1,3
+ do i=1,n
+ ind=(i-1)*3+k
+ d_t_work(ind)=0.0d0
+ do j=1,n
+ d_t_work(ind)=d_t_work(ind)
+ & +Gvec(i,j)*d_t_work_new((j-1)*3+k)
+ enddo
+c write (iout,*) "i",i," ind",ind," d_t_work",d_t_work(ind)
+c call flush(iout)
+ enddo
+ enddo
+#ifdef DEBUG
+ aux=0.0d0
+ do k=1,3
+ do i=1,n
+ do j=1,n
+ aux=aux+inertia(i,j)*d_t_work(3*(i-1)+k)*d_t_work(3*(j-1)+k)
+ enddo
+ enddo
+ enddo
+ Ek3=Ek3+aux/2
+#endif
+c Transfer to the d_t vector
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+ ind=0
+c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+ do i=innt,inct
+ do j=1,3
+ ind=ind+1
+ d_t(j,i)=d_t_work(ind)
+ enddo
+ if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
+ do j=1,3
+ ind=ind+1
+ d_t(j,i+nres)=d_t_work(ind)
+ enddo
+ endif
+ enddo
+ enddo
+ if (large) then
+ write (iout,*)
+ write (iout,*) "Random velocities in the Calpha,SC space"
+ do i=1,nres
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
+ & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ endif
+ call kinetic_CASC(Ek1)
+!
+! Transform the velocities to virtual-bond space
+!
+#define WLOS
+#ifdef WLOS
+ if (nnt.eq.1) then
+ d_t(:,0)=d_t(:,1)
+ endif
+ do i=1,nres
+ if (itype(i).eq.10 .or. itype(i).eq.ntyp1) then
+ do j=1,3
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ else
+ do j=1,3
+ d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ end if
+ enddo
+ d_t(:,nres)=0.0d0
+ d_t(:,nct)=0.0d0
+ d_t(:,2*nres)=0.0d0
+ if (nnt.gt.1) then
+ d_t(:,0)=d_t(:,1)
+ d_t(:,1)=0.0d0
+ endif
+c d_a(:,0)=d_a(:,1)
+c d_a(:,1)=0.0d0
+c write (iout,*) "Shifting accelerations"
+ do ichain=2,nchain
+c write (iout,*) "ichain",chain_border1(1,ichain)-1,
+c & chain_border1(1,ichain)
+ d_t(:,chain_border1(1,ichain)-1)=d_t(:,chain_border1(1,ichain))
+ d_t(:,chain_border1(1,ichain))=0.0d0
+ enddo
+c write (iout,*) "Adding accelerations"
+ do ichain=2,nchain
+c write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+c & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=
+ & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+ do ichain=2,nchain
+ write (iout,*) "chain",ichain,chain_border1(1,ichain)-1,
+ & chain_border(2,ichain-1)
+ d_t(:,chain_border1(1,ichain)-1)=
+ & d_t(:,chain_border1(1,ichain)-1)+d_t(:,chain_border(2,ichain-1))
+ d_t(:,chain_border(2,ichain-1))=0.0d0
+ enddo
+#else
+ ibeg=0
+c do j=1,3
+c d_t(j,0)=d_t(j,nnt)
+c enddo
+ do ichain=1,nchain
+ innt=chain_border(1,ichain)
+ inct=chain_border(2,ichain)
+c write (iout,*) "ichain",ichain," innt",innt," inct",inct
+c write (iout,*) "ibeg",ibeg
+ do j=1,3
+ d_t(j,ibeg)=d_t(j,innt)
+ enddo
+ ibeg=inct+1
+ do i=innt,inct
+ if (iabs(itype(i).eq.10)) then
+c write (iout,*) "i",i,(d_t(j,i),j=1,3),(d_t(j,i+1),j=1,3)
+ do j=1,3
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ else
+ do j=1,3
+ d_t(j,i+nres)=d_t(j,i+nres)-d_t(j,i)
+ d_t(j,i)=d_t(j,i+1)-d_t(j,i)
+ enddo
+ end if
+ enddo
+ enddo
+#endif
+ if (large) then
+ write (iout,*)
+ write (iout,*)
+ & "Random velocities in the virtual-bond-vector space"
+ write (iout,'(3hORG,1h(,i5,1h),3f10.5)') 0,(d_t(j,0),j=1,3)
+ do i=1,nres
+ write (iout,'(a3,1h(,i5,1h),3f10.5,3x,3f10.5)')
+ & restyp(itype(i)),i,(d_t(j,i),j=1,3),(d_t(j,i+nres),j=1,3)
+ enddo
+ write (iout,*)
+ write (iout,*) "Kinetic energy from inertia matrix eigenvalues",
+ & Ek
+ write (iout,*)
+ & "Kinetic temperatures from inertia matrix eigenvalues",
+ & 2*Ek/(3*dimen*Rb)
+#ifdef DEBUG
+ write (iout,*) "Kinetic energy from inertia matrix",Ek3
+ write (iout,*) "Kinetic temperatures from inertia",
+ & 2*Ek3/(3*dimen*Rb)
+#endif
+ write (iout,*) "Kinetic energy from velocities in CA-SC space",
+ & Ek1
+ write (iout,*)
+ & "Kinetic temperatures from velovities in CA-SC space",
+ & 2*Ek1/(3*dimen*Rb)
+ call kinetic(Ek1)
+ write (iout,*)
+ & "Kinetic energy from virtual-bond-vector velocities",Ek1
+ write (iout,*)
+ & "Kinetic temperature from virtual-bond-vector velocities ",
+ & 2*Ek1/(dimen3*Rb)
+ endif
+#else