********>Amber 10 Update for trajectory post processing Author: Daniel R. Roe, 06-18-2009 Programs: AMBER 10 Description: Modify post processing of trajectory files (imin=5). Will now work for unformatted input trajectories (ioutfm=1) and PBC (ntb!=0). This will also update the trajene test case and change the sander/Makefile to force a rebuild of the depend file on a 'make clean'. Requires that all previous bugfixes (1-13) are applied and a 'make clean' in order to compile. Usage: Save this file in your $AMBERHOME directory and then apply this patch file to your amber 10 distribution as follows: cd $AMBERHOME patch -p0 -N -r patch_rejects < amber10.imin5.patch ----------------------------------------------------------------------------- --- src/sander/trajene.f 2009-06-17 16:10:58.000000000 -0400 +++ src/sander/trajene.f 2009-06-16 17:29:05.000000000 -0400 @@ -1,6 +1,19 @@ +! #include "copyright.h" #include "dprec.h" #include "assert.h" + +module trajenemod + +! MODULE: TRAJENEMOD +! ================== TRAJECTORY ENERGY POST-PROCESSING =================== +! Daniel R. Roe, 2009 +! Based on original implementation by Carlos Simmerling + +implicit none + +contains + !********************************************************************* ! SUBROUTINE TRAJENE !********************************************************************* @@ -9,51 +22,114 @@ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine trajene here] subroutine trajene(x,ix,ih,ipairs,ene,ok,qsetup) - - use nblist,only: fill_tranvec,a,b,c!,alpha,beta,gamma,nbflag,cutoffnb + !\/-----\/--For Debug + use nblist,only: fill_tranvec,a,b,c,alpha,beta,gamma!,nbflag,cutoffnb use constants,only: one,half +#ifdef BINTRAJ + use netcdf + use bintraj,only: end_binary_frame +#endif implicit none - integer ipairs(*) - _REAL_ carrms - - _REAL_ x(*),ene(*),oldbox(3),rmu(3) - integer ix(*) + ! INPUT VARIABLES + integer ipairs(*),ix(*) + _REAL_ x(*),ene(*) character(len=4) ih(*) - logical ok,qsetup,loutfm - integer member,j + logical ok,qsetup + + ! INTERNAL VARIABLES + integer member,j,xstop + integer newnfft1,newnfft2,newnfft3 + _REAL_ carrms,oldbox(3),rmu(3) + logical loutfm + +#ifdef BINTRAJ + ! For netcdf files + integer ncid,frameID,ncframe,atomID,ncatom,coordID,cellLengthID,cellAngleID + integer err +#endif # include "memory.h" # include "tgtmd.h" ! DAN ROE: Is extra.h still needed? # include "extra.h" ! needed for ntb, box -# include "box.h" +# include "box.h" +! Needed for debug info !# include "ew_cntrl.h" !# include "ew_erfc_spline.h" ! needed for nfft # include "ew_pme_recip.h" ! Needed for INPTRAJ_UNIT, inptraj # include "files.h" - +!------------------------------------------------ loutfm = ioutfm <= 0 - ! Open the inptraj file for coordinate reading - call amopen(INPTRAJ_UNIT,inptraj,'O','F','R') - - read(INPTRAJ_UNIT,100) title - write (6,100) title - 100 format(a80) + ! If PBC, save original nfft sizes. PME allocates memory for nfft based + ! on the original box size, so if the box becomes much larger than the + ! original box size the new nfft values would cause a memory error. + if (ntb>0) then + write (6,'(a,3i6)') "TRAJENE: Original NFFTs: ",nfft1,nfft2,nfft3 + endif - member=0 + ! This is for indexing the x array to get coordinates + xstop=lcrd+(natom*3)-1; - ! loop over trajectory file, exiting only on error or end of file + ! Open the inptraj file for coordinate reading + if (ioutfm==1) then +#ifdef BINTRAJ + ! Open NETCDF format trajectory file + err = nf90_open(inptraj,nf90_nowrite,ncid) + if (err /= nf90_noerr) then + write(6,'(a,a)') "TRAJENE: Could not open netcdf file: ",inptraj + return + endif + ! Get netcdf file information + ! NOTE: NEEDS BETTER ERROR HANDLING + err = nf90_inq_dimid(ncid,"frame",frameID) + err = nf90_inquire_dimension(ncid,frameID,len = ncframe) + write (6,'(a,i)') "TRAJENE: Frames in trajectory= ",ncframe + err = nf90_inq_dimid(ncid,"atom",atomID) + err = nf90_inquire_dimension(ncid,atomID,len = ncatom) + write (6,'(a,i)') "TRAJENE: Atoms in trajectory= ",ncatom + if (ncatom /= natom) then + write (6,'(a)') "TRAJENE: # atoms in traj does not match # atoms in top." + return + endif + err = nf90_inq_varid(ncid,"coordinates",coordID) + if (ntb>0) then + err = nf90_inq_varid(ncid,"cell_lengths",cellLengthID) + err = nf90_inq_varid(ncid,"cell_angles",cellAngleID) + endif +#else + ! No NETCDF support + write(6,'(a)') "TRAJENE: Netcdf support not compiled in this version. Please recompile with BINTRAJ" + return +#endif + else + ! Standard Amber Trajectory Format + call amopen(INPTRAJ_UNIT,inptraj,'O','F','R') + read(INPTRAJ_UNIT,'(a80)') title + write (6,'(a80)') title + endif + ! Now loop over trajectory file, exiting only on error or end of file + member=1 do while ( .true. ) ! --- read next coordinate set from trajectory - - read(INPTRAJ_UNIT,110,end=1000,err=1010) (x(j),j=lcrd,lcrd+natom*3-1) + if (ioutfm/=1) then + read(INPTRAJ_UNIT,'(10f8.3)',end=1000,err=1010) (x(j),j=lcrd,xstop) +#ifdef BINTRAJ + else + if (member>ncframe) goto 1000 + + err = nf90_get_var(ncid,coordID,x(lcrd:xstop), & + start = (/ 1, 1, member /), & + count = (/ 3, (xstop-lcrd+1)/3, 1 /) ) + if (err /= nf90_noerr) goto 1010 +#endif + endif ! DAN ROE: ! If box coords are present in trajectory (ntb>0), read them in and @@ -67,33 +143,43 @@ !write(6,*) "DEBUG: OLDabc: ",a,b,c ! 2- Read in current box coords. - read(INPTRAJ_UNIT,'(3f8.3)',end=1000,err=1010) box(1), box(2), box(3) + if (ioutfm/=1) then + read(INPTRAJ_UNIT,'(3f8.3)',end=1000,err=1020) box(1), box(2), box(3) +#ifdef BINTRAJ + else + err = nf90_get_var(ncid,cellLengthID,box(1:3), & + start = (/ 1, member /), & + count = (/ 3, 1 /) ) + if (err /= nf90_noerr) goto 1020 +#endif + endif !write(6,*) "DEBUG: BOX: ",box(1),box(2),box(3) ! 3- If the box size has changed, some parameters have to be recalc. if (box(1)/=oldbox(1).or.box(2)/=oldbox(2).or.box(3)/=oldbox(3)) then - ! 3a- Calculate box scaling factors - rmu(1) = box(1) / oldbox(1) - rmu(2) = box(2) / oldbox(2) - rmu(3) = box(3) / oldbox(3) - !write(6,*) "DEBUG: RMU: ",rmu(1),rmu(2),rmu(3) - - ! 3b- Update the box - ! DAN ROE: Make sure fill_tranvec is needed. - call redo_ucell(rmu) - call fill_tranvec() - !write(6,*) "DEBUG: NEWabc: ",a,b,c + !write (6,'(a)') "TRAJENE: Updating Box parameters" + ! 3.a - Update PME cell parameters + call fill_ucell(box(1),box(2),box(3),alpha,beta,gamma) + !call fill_tranvec() ! 3c- Recompute grid sizes + ! If the grid sizes change the non-bonded energy will probably be + ! significantly different than it was in the simulation. + ! ! RCFFT needs an even dimension for x direction - call compute_nfft((a + one)*half ,nfft1) - nfft1=nfft1*2 - call compute_nfft(b,nfft2) - call compute_nfft(c,nfft3) + call compute_nfft((a + one)*half ,newnfft1) + newnfft1=newnfft1*2 + call compute_nfft(b,newnfft2) + call compute_nfft(c,newnfft3) + !write (6,'(a,3i6)') "TRAJENE: New NFFTs: ",newnfft1,newnfft2,newnfft3 + if ( (newnfft1/=nfft1) .or. (newnfft2/=nfft2) .or. (newnfft3/=nfft3) ) then + write(6,'(a)') 'TRAJENE WARNING: nfft size would change with this box!' + write (6,'(a,3i6)') "TRAJENE: New NFFTs: ",newnfft1,newnfft2,newnfft3 + endif endif ! box sizes have changed - ! DAN ROE: More Debug + ! DAN ROE: More Debug for PME !write(6,'(/a)') 'Ewald parameters:' !write(6,'(5x,4(a,i8))') 'verbose =',verbose, & ! ', ew_type =',ew_type,', nbflag =',nbflag, & @@ -114,11 +200,8 @@ !9007 format (5x,'Ewald Coefficient =',f9.5) endif ! ntb>0 - 110 format(10f8.3) - - member=member+1 - write (6,'(a,i6)') 'minimizing coord set #',member + member=member+1 call runmin(x,ix,ih,ipairs,x(lcrd),x(lforce),x(lvel), & ix(iibh),ix(ijbh),x(l50),x(lwinv),ix(ibellygp), & @@ -134,10 +217,14 @@ ! write the frame to the mdcrd file if the user has set ntwx=1 ! don't worry about imaging etc since we write the same way it came in + ! NOTE: Eventually make ntwx work properly, i.e. frames can be skipped if (master .and. ntwx >= 1) then call corpac(x(lcrd),1,natom*3,MDCRD_UNIT,loutfm) - if(ntb > 0) call corpac(box,1,3,MDCRD_UNIT,loutfm) + if (ntb > 0) call corpac(box,1,3,MDCRD_UNIT,loutfm) +#ifdef BINTRAJ + if (ioutfm == 1) call end_binary_frame(MDCRD_UNIT) +#endif !elseif (master) then ! write (6,*) "Not writing coordinates to mdcrd due to NTWX value" endif @@ -148,13 +235,22 @@ ! ---end of trajectory file - 1000 write (6,*) "Trajectory file ended" + 1000 write (6,'(a)') "TRAJENE: Trajectory file ended" ok=.true. - return + goto 1500 + + 1010 write (6,'(a)') "TRAJENE: Error reading trajectory coordinates" + ok=.false. + goto 1500 - 1010 write (6,*) "Error in trajectory file" + 1020 write (6,'(a)') "TRAJENE: Error reading box coordinates" ok=.false. + goto 1500 + + 1500 write (6,'(a)') "TRAJENE: Trajene complete." return end subroutine trajene + +end module trajenemod --- src/sander/sander.f 2009-06-17 16:10:55.000000000 -0400 +++ src/sander/sander.f 2009-06-17 16:19:50.000000000 -0400 @@ -71,6 +71,9 @@ #endif /* MPI */ use pimd_vars, only: ipimd,ineb + + use trajenemod, only: trajene + !!jtc ========================= PUPIL INTERFACE ========================= #ifdef PUPIL_SUPPORT ! Using data structure for PUPIL --- src/sander/mdread.f 2009-06-17 16:10:54.000000000 -0400 +++ src/sander/mdread.f 2009-06-17 16:21:05.000000000 -0400 @@ -1729,13 +1729,11 @@ write(6,'(/2x,a,i3,a)') 'IMIN (',imin,') must be >= 0.' inerr = 1 end if - if (imin == 5 .and. ioutfm /= 0 ) then - write(6,'(/2x,a,i3,a)') 'IMIN=5 currently requires formatted mdcrd/inptraj (ioutfm=0).' - inerr = 1 - end if - if (imin == 5 .and. ifbox /= 0 .and. ntb /= 0) then - write(6,'(/2x,a,i3,a)') 'IMIN=5 does not support periodic boundaries (ifbox>0, ntb>0).' - inerr = 1 + if (imin == 5) then + if (ifbox /= 0 .and. ntb /= 0) then + write(6,'(/2x,a)') 'WARNING: IMIN=5 with periodic boundaries can result in odd energies being calculated' + write(6,'(/2x,a)') ' when box size can vary (i.e. ntb=2). Use with caution.' + endif end if --- test/trajene/trajene.out.save 2009-06-18 09:38:48.000000000 -0400 +++ test/trajene/trajene.out.save 2009-06-18 09:50:48.000000000 -0400 @@ -2156,9 +2156,10 @@ Current target RMSD: 0.000 minimization completed, ENE= -.29125793E+03 RMS= 0.125516E+02 Final RMSD from reference: 2.539 - Trajectory file ended +TRAJENE: Trajectory file ended +TRAJENE: Trajene complete. -------------------------------------------------------------------------------- 5. TIMINGS -------------------------------------------------------------------------------- --- src/sander/Makefile 2009-06-18 10:56:31.000000000 -0400 +++ src/sander/Makefile 2009-06-18 10:57:32.000000000 -0400 @@ -283,6 +283,7 @@ -/bin/rm *.mod -/bin/rm *.d -/bin/rm *nbflag + /bin/rm depend; make depend cd ../dcqtp; make clean depend::