! #include "copyright.h" #include "../include/dprec.fh" #include "../include/assert.fh" module trajenemod ! MODULE: TRAJENEMOD ! ================== TRAJECTORY ENERGY POST-PROCESSING =================== ! Daniel R. Roe, 2009 ! Based on original implementation by Carlos Simmerling implicit none contains !********************************************************************* ! SUBROUTINE TRAJENE !********************************************************************* ! carlos add trajene routine for processing trajectory energies !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+ [Enter a one-line description of subroutine trajene here] subroutine trajene(x,ix,ih,ipairs,ene,ok,qsetup) !\/-----\/--For Debug use nblist,only: fill_tranvec,a,b,c,alpha,beta,gamma!,nbflag,cutoffnb use constants,only: one,half use state use file_io_dat, only : INPTRAJ_UNIT, inptraj, ioutfm, MDCRD_UNIT, & ntwx, title #ifdef BINTRAJ use netcdf use bintraj,only: end_binary_frame #endif implicit none ! INPUT VARIABLES integer ipairs(*),ix(*) _REAL_ x(*) type(state_rec) :: ene character(len=4) ih(*) 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" ! Needed for debug info !# include "ew_cntrl.h" !# include "ew_erfc_spline.h" ! needed for nfft # include "ew_pme_recip.h" !------------------------------------------------ loutfm = ioutfm <= 0 ! 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 ! This is for indexing the x array to get coordinates xstop=lcrd+(natom*3)-1; ! 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 write(6,'(a,a)') "TRAJENE: ",trim(nf90_strerror(err)) 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,i5)') "TRAJENE: Frames in trajectory= ",ncframe err = nf90_inq_dimid(ncid,"atom",atomID) err = nf90_inquire_dimension(ncid,atomID,len = ncatom) write (6,'(a,i5)') "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 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 ! update the unit cell and grid size. if (ntb>0) then ! 1- Save old box coords. oldbox(1) = box(1) oldbox(2) = box(2) oldbox(3) = box(3) !write(6,*) "DEBUG: OLDBOX: ",oldbox(1),oldbox(2),oldbox(3) !write(6,*) "DEBUG: OLDabc: ",a,b,c ! 2- Read in current box coords. 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 !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 ,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 for PME !write(6,'(/a)') 'Ewald parameters:' !write(6,'(5x,4(a,i8))') 'verbose =',verbose, & ! ', ew_type =',ew_type,', nbflag =',nbflag, & ! ', use_pme =',use_pme !write(6,'(5x,4(a,i8))') 'vdwmeth =',vdwmeth, & ! ', eedmeth =',eedmeth,', netfrc =',netfrc !write(6, 9002) a, b, c !write(6, 9003) alpha, beta, gamma !write(6, 9004) nfft1, nfft2, nfft3 !write(6, 9006) cutoffnb, dsum_tol !write(6, 9007) ew_coeff !write(6, 9005) order !9002 format (5x,'Box X =',f9.3,3x,'Box Y =',f9.3,3x,'Box Z =',f9.3) !9003 format (5x,'Alpha =',f9.3,3x,'Beta =',f9.3,3x,'Gamma =',f9.3) !9004 format (5x,'NFFT1 =',i5 ,7x,'NFFT2 =',i5 ,7x,'NFFT3 =',i5) !9005 format (5x,'Interpolation order =',i5) !9006 format (5x,'Cutoff=',f9.3,3x,'Tol =',e9.3) !9007 format (5x,'Ewald Coefficient =',f9.5) endif ! ntb>0 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), & x(l95),ene,carrms,qsetup) write (6,364) ene%pot%tot,carrms 364 format ('minimization completed, ENE=',1x,e14.8, & 1x,'RMS=',1x,e12.6) if (master .and. itgtmd == 1) then write (6,'(a,f8.3)') "Final RMSD from reference: ",rmsdvalue end if ! 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) #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 ! ---loop for next coordinate set end do ! ---end of trajectory file 1000 write (6,'(a)') "TRAJENE: Trajectory file ended" ok=.true. goto 1500 1010 write (6,'(a)') "TRAJENE: Error reading trajectory coordinates" ok=.false. goto 1500 1020 write (6,'(a)') "TRAJENE: Error reading box coordinates" ok=.false. goto 1500 1500 write (6,'(a)') "TRAJENE: Trajene complete." #ifdef BINTRAJ if (ioutfm==1) then err=nf90_close(ncid) if (err /= nf90_noerr) then write(6,'(a,a)') "TRAJENE: Could not close netcdf file: ",inptraj write(6,'(a,a)') "TRAJENE: ",trim(nf90_strerror(err)) endif endif #endif return end subroutine trajene end module trajenemod