#ifdef FIVEDIAG subroutine kinetic(KE_total) c---------------------------------------------------------------- c This subroutine calculates the total kinetic energy of the chain c----------------------------------------------------------------- c 3/5/2020 AL Corrected for multichain systems, no fake peptide groups c inside, implemented with five-diagonal inertia matrix implicit none include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' include 'COMMON.LAGRANGE.5diag' include 'COMMON.IOUNITS' double precision KE_total integer i,j,k,iti double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3), & mag1,mag2,v(3) KEt_p=0.0d0 KEt_sc=0.0d0 KEr_p=0.0D0 KEr_sc=0.0D0 c write (iout,*) "ISC",(isc(itype(i)),i=1,nres) c The translational part for peptide virtual bonds do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3 c Skip dummy peptide groups if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo c write (iout,*) "Kinetic trp:",i,(v(j),j=1,3) vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) endif do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo c write(iout,*) 'KEt_p', KEt_p c The translational part for the side chain virtual bond c Only now we can initialize incr with zeros. It must be equal c to the velocities of the first Calpha. do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct iti=iabs(itype(i)) if (itype(i).eq.10 .and. itype(i).ne.ntyp1) then do j=1,3 v(j)=incr(j) enddo else do j=1,3 v(j)=incr(j)+d_t(j,nres+i) enddo endif c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) c 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)) 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) enddo enddo c goto 111 c write(iout,*) 'KEt_sc', KEt_sc c The part due to stretching and rotation of the peptide groups do i=nnt,nct-1 if (itype(i).ne.ntyp1.and.itype(i+1).ne.ntyp1) then c write (iout,*) "i",i c write (iout,*) "i",i," mag1",mag1," mag2",mag2 do j=1,3 incr(j)=d_t(j,i) enddo c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) & +incr(3)*incr(3)) endif enddo c goto 111 c write(iout,*) 'KEr_p', KEr_p c The rotational part of the side chain virtual bond do i=nnt,nct iti=iabs(itype(i)) if (itype(i).ne.10.and.itype(i).ne.ntyp1) then do j=1,3 incr(j)=d_t(j,nres+i) enddo c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ & incr(3)*incr(3)) endif enddo c The total kinetic energy 111 continue c write(iout,*) ' KEt_p',KEt_p,' KEt_sc',KEt_sc,' KEr_p',KEr_p, c & ' KEr_sc', KEr_sc KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc) c write (iout,*) "KE_total",KE_total return end #else subroutine kinetic(KE_total) c---------------------------------------------------------------- c This subroutine calculates the total kinetic energy of the chain c----------------------------------------------------------------- implicit none include 'DIMENSIONS' include 'COMMON.VAR' include 'COMMON.CHAIN' include 'COMMON.DERIV' include 'COMMON.GEO' include 'COMMON.LOCAL' include 'COMMON.INTERACT' include 'COMMON.MD' include 'COMMON.LAGRANGE' include 'COMMON.IOUNITS' double precision KE_total integer i,j,k,iti double precision KEt_p,KEt_sc,KEr_p,KEr_sc,incr(3), & mag1,mag2,v(3) KEt_p=0.0d0 KEt_sc=0.0d0 c write (iout,*) "ISC",(isc(itype(i)),i=1,nres) c The translational part for peptide virtual bonds do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct-1 c write (iout,*) "Kinetic trp:",i,(incr(j),j=1,3) do j=1,3 v(j)=incr(j)+0.5d0*d_t(j,i) enddo vtot(i)=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) KEt_p=KEt_p+(v(1)*v(1)+v(2)*v(2)+v(3)*v(3)) do j=1,3 incr(j)=incr(j)+d_t(j,i) enddo enddo c write(iout,*) 'KEt_p', KEt_p c The translational part for the side chain virtual bond c Only now we can initialize incr with zeros. It must be equal c to the velocities of the first Calpha. do j=1,3 incr(j)=d_t(j,0) enddo do i=nnt,nct iti=iabs(itype(i)) if (itype(i).eq.10) then do j=1,3 v(j)=incr(j) enddo else do j=1,3 v(j)=incr(j)+d_t(j,nres+i) enddo endif c write (iout,*) "Kinetic trsc:",i,(incr(j),j=1,3) c 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)) 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) enddo enddo c goto 111 c write(iout,*) 'KEt_sc', KEt_sc c The part due to stretching and rotation of the peptide groups KEr_p=0.0D0 do i=nnt,nct-1 c write (iout,*) "i",i c write (iout,*) "i",i," mag1",mag1," mag2",mag2 do j=1,3 incr(j)=d_t(j,i) enddo c write (iout,*) "Kinetic rotp:",i,(incr(j),j=1,3) KEr_p=KEr_p+(incr(1)*incr(1)+incr(2)*incr(2) & +incr(3)*incr(3)) enddo c goto 111 c write(iout,*) 'KEr_p', KEr_p c The rotational part of the side chain virtual bond KEr_sc=0.0D0 do i=nnt,nct iti=iabs(itype(i)) if (itype(i).ne.10) then do j=1,3 incr(j)=d_t(j,nres+i) enddo c write (iout,*) "Kinetic rotsc:",i,(incr(j),j=1,3) KEr_sc=KEr_sc+Isc(iti)*(incr(1)*incr(1)+incr(2)*incr(2)+ & incr(3)*incr(3)) endif enddo c The total kinetic energy 111 continue c write(iout,*) 'KEr_sc', KEr_sc KE_total=0.5d0*(mp*KEt_p+KEt_sc+0.25d0*Ip*KEr_p+KEr_sc) c write (iout,*) "KE_total",KE_total return end #endif