********> bugfix.78 Author: Gerhard Hummer (Los Alamos) and Bill Ross Date: 11/2/96 Programs: Sander Severity: Moderate Problem: SHAKE failure in Ewald runs. Affects: Ewald sum calculations. Cause: In Ewald runs, the box coordinates read from mdin (NTX<7) or restrt (NTX=7) files are not copied to the general box variables until the first printing of coords to mdcrd (step=NTWX). This results in prmtop box coords being used for the general box variables, which are used when calculating box-crossing, bond-related terms, while the Ewald box variables are used for adjusting the box when pressure scaling. Thus if these are out of sync, SHAKE failure can occur, and in general errors would creep into the calculation. Also, when the box is scaled in constant pressure Ewald, the general box is updated. Fix: Make the following changes: ---------------------------------------------------------------------------- Replace box.h with the following: ---------------------------------------------------------------------------- c c ... floats: c #ifdef DPREC double precision BOX,BOXI,BETA,BOXH,BOXHM,BOXHM2, + BOXOH,BOXOQ,COSB,COSB2,BOCT,fothi, + CUT,CUT2ND,SCNB,SCEE,DIELC, + CUTCAP,XCAP,YCAP,ZCAP,FCAP, + RWELL,RAD,WEL,RADHB,WELHB #else real BOX,BOXI,BETA,BOXH,BOXHM,BOXHM2, + BOXOH,BOXOQ,COSB,COSB2,BOCT,fothi, + CUT,CUT2ND,SCNB,SCEE,DIELC, + CUTCAP,XCAP,YCAP,ZCAP,FCAP, + RWELL,RAD,WEL,RADHB,WELHB #endif COMMON/BOXR/BOX(3),BOXI(3),BETA,BOXH(3),BOXHM,BOXHM2, + BOXOH,BOXOQ,COSB,COSB2,BOCT(3),fothi, + CUT,CUT2ND,SCNB,SCEE,DIELC, + CUTCAP,XCAP,YCAP,ZCAP,FCAP, + RWELL,RAD(100),WEL(100),RADHB(200),WELHB(200) c c ... integers: c integer NTM,NTB,ifbox,IDIEL,NBUCK,NUMPK,NBIT, + IFTRES,IPTRES,IPTATM,NSPSOL,NSPSTR,IPTSOL,IFCRST, + IMGSLT,nbsola,nbsolh,nsolu,nsolv, + IFCAP,NATCAP,ISFTRP,ISFT01,jwflag c COMMON/BOXI/NTM,NTB,ifbox,IDIEL,NBUCK,NUMPK,NBIT, + IFTRES,IPTRES,IPTATM,NSPSOL,NSPSTR,IPTSOL,IFCRST, + IMGSLT,nbsola,nbsolh,nsolu,nsolv, + IFCAP,NATCAP,ISFTRP,ISFT01,jwflag ---------------------------------------------------------------------------- Replace ewald/unitcell.h with the following: ---------------------------------------------------------------------------- c c DEFINING QUANTITIES FOR UNIT CELL c ucell is the 3x3 of direct lattice vectors c dirlng are their lengths c recip is the 3x3 of reciprocal lattice vectors c reclng are their lengths c--------------------------------------------------------- _REAL_ a,b,c,alpha,betta,gamma,volume _REAL_ ucell,recip,dirlng,reclng,sphere,cutoffnb _REAL_ olducell,oldrecip common/unitcell/ucell(3,3),recip(3,3),dirlng(3),reclng(3), $ olducell(3,3),oldrecip(3,3), $ a,b,c,alpha,betta,gamma,volume,sphere,cutoffnb ---------------------------------------------------------------------------- Make the following changes to getcor.f: ---------------------------------------------------------------------------- *** OLD getcor.f --- NEW getcor.f *************** *** 114,119 **** --- 114,127 ---- endif call fill_ucell(a,b,c,alpha,betta,gamma) write(6,9129) a,b,c,alpha,betta,gamma + c + c -- set ewald-independent versions of box variables + c + box(1) = a + box(2) = b + box(3) = c + beta = betta + c else c c -- normal case, not running ewald ---------------------------------------------------------------------------- Make the following changes to mdread.f: ---------------------------------------------------------------------------- *** OLD mdread.f --- NEW mdread.f *************** *** 551,557 **** c-------------------------------------------------------------------- c if ( iewald .eq. 1 ) then ! call load_ewald_info(5,cut) endif c c =========================== END EWALD =========================== --- 551,557 ---- c-------------------------------------------------------------------- c if ( iewald .eq. 1 ) then ! call load_ewald_info(5) endif c c =========================== END EWALD =========================== ---------------------------------------------------------------------------- Make the following changes to ewald.f: ---------------------------------------------------------------------------- *** OLD ewald.f --- NEW ewald.f *************** *** 1902,1918 **** c------------------------------------------------------------------- c --- FILL_UCELL --- ! subroutine fill_ucell(ta,tb,tc,talpha,tbeta,tgamma) implicit none # include "ewald/unitcell.h" ! _REAL_ ta,tb,tc,talpha,tbeta,tgamma a = ta b = tb c = tc alpha = talpha ! beta = tbeta gamma = tgamma ! call get_ucell(a,b,c,alpha,beta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume) return end --- 1902,1922 ---- c------------------------------------------------------------------- c --- FILL_UCELL --- ! subroutine fill_ucell(ta,tb,tc,talpha,tbetta,tgamma) implicit none # include "ewald/unitcell.h" ! # include "box.h" ! _REAL_ ta,tb,tc,talpha,tbetta,tgamma a = ta b = tb c = tc + box(1) = a + box(2) = b + box(3) = c alpha = talpha ! betta = tbetta gamma = tgamma ! call get_ucell(a,b,c,alpha,betta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume) return end *************** *** 1972,1986 **** c ...simple wrapper to set the values of the arguments to c the unit cell parameters... c ! subroutine find_ucell(ta,tb,tc,talpha,tbeta,tgamma) implicit none # include "ewald/unitcell.h" ! _REAL_ ta,tb,tc,talpha,tbeta,tgamma ta = a tb = b tc = c talpha = alpha ! tbeta = beta tgamma = gamma return end --- 1976,1990 ---- c ...simple wrapper to set the values of the arguments to c the unit cell parameters... c ! subroutine find_ucell(ta,tb,tc,talpha,tbetta,tgamma) implicit none # include "ewald/unitcell.h" ! _REAL_ ta,tb,tc,talpha,tbetta,tgamma ta = a tb = b tc = c talpha = alpha ! tbetta = betta tgamma = gamma return end *************** *** 2408,2414 **** c in the x-y plane with positive y, and that the direct lattice c vectors are a non-degenerate right handed system. Thus the 3rd c vector has positive z component. Alpha is the angle (in degrees) ! c between 2nd and 3rd vectors, beta is the angle (in degrees) c between 1st and 3rd vectors, and gamma is the angle (in degrees) c between 1st and 2nd vectors. The lengths of the direct lattice c vectors are given by dirlng(1), dirlng(2) and dirlng(3), --- 2412,2418 ---- c in the x-y plane with positive y, and that the direct lattice c vectors are a non-degenerate right handed system. Thus the 3rd c vector has positive z component. Alpha is the angle (in degrees) ! c between 2nd and 3rd vectors, betta is the angle (in degrees) c between 1st and 3rd vectors, and gamma is the angle (in degrees) c between 1st and 2nd vectors. The lengths of the direct lattice c vectors are given by dirlng(1), dirlng(2) and dirlng(3), *************** *** 2415,2424 **** c whereas the reciprocal lengths of the reciprocal vectors are c reclng(1),reclng(2) and reclng(3). c ! subroutine get_ucell(a,b,c,alpha,beta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume) implicit none ! _REAL_ a,b,c,alpha,beta,gamma _REAL_ sphere,volume _REAL_ ucell(3,3),recip(3,3),dirlng(3),reclng(3) _REAL_ factor,u23(3),u31(3),u12(3) --- 2419,2428 ---- c whereas the reciprocal lengths of the reciprocal vectors are c reclng(1),reclng(2) and reclng(3). c ! subroutine get_ucell(a,b,c,alpha,betta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume) implicit none ! _REAL_ a,b,c,alpha,betta,gamma _REAL_ sphere,volume _REAL_ ucell(3,3),recip(3,3),dirlng(3),reclng(3) _REAL_ factor,u23(3),u31(3),u12(3) *************** *** 2432,2438 **** ucell(1,2) = b*cos(factor*gamma) ucell(2,2) = b*sin(factor*gamma) ucell(3,2) = 0.d0 ! ucell(1,3) = c*cos(factor*beta) ucell(2,3) = $ (b*c*cos(factor*alpha)-ucell(1,3)*ucell(1,2))/ucell(2,2) ucell(3,3) = sqrt(c**2-ucell(1,3)**2-ucell(2,3)**2) --- 2436,2442 ---- ucell(1,2) = b*cos(factor*gamma) ucell(2,2) = b*sin(factor*gamma) ucell(3,2) = 0.d0 ! ucell(1,3) = c*cos(factor*betta) ucell(2,3) = $ (b*c*cos(factor*alpha)-ucell(1,3)*ucell(1,2))/ucell(2,2) ucell(3,3) = sqrt(c**2-ucell(1,3)**2-ucell(2,3)**2) *************** *** 2940,2951 **** c for locmem() and must be called prior to calling locmem(). c Currently this is called inside mdread(). c ! subroutine load_ewald_info(nf,cut) implicit none # include "ewald/ewcntrl.h" # include "ewald/unitcell.h" # include "ewald/pme_recip.h" ! _REAL_ cut integer nf write(6, 9001) --- 2944,2955 ---- c for locmem() and must be called prior to calling locmem(). c Currently this is called inside mdread(). c ! subroutine load_ewald_info(nf) implicit none # include "ewald/ewcntrl.h" # include "ewald/unitcell.h" # include "ewald/pme_recip.h" ! # include "box.h" integer nf write(6, 9001) *************** *** 2967,2977 **** #endif call start_numtasks() cutoffnb = cut ! call read_ewald(nf,a,b,c,alpha,beta,gamma, $ nfft1,nfft2,nfft3,order,ischrgd,verbose, $ checkacc,dsum_tol) 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, 9005) order --- 2971,2985 ---- #endif call start_numtasks() cutoffnb = cut ! call read_ewald(nf,a,b,c,alpha,betta,gamma, $ nfft1,nfft2,nfft3,order,ischrgd,verbose, $ checkacc,dsum_tol) + box(1) = a + box(2) = b + box(3) = c + beta = betta write(6, 9002) a, b, c ! write(6, 9003) alpha, betta, gamma write(6, 9004) nfft1, nfft2, nfft3 write(6, 9006) cutoffnb, dsum_tol write(6, 9005) order *************** *** 2986,2992 **** endif c c get some related quantities ! call get_ucell(a,b,c,alpha,beta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume) call find_ewaldcof(cutoffnb,dsum_tol,ew_coeff) write(6, 9007) ew_coeff --- 2994,3000 ---- endif c c get some related quantities ! call get_ucell(a,b,c,alpha,betta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume) call find_ewaldcof(cutoffnb,dsum_tol,ew_coeff) write(6, 9007) ew_coeff *************** *** 3674,3690 **** c------------------------------------------------------------------- c --- READ_EWALD --- ! subroutine read_ewald(nf,a,b,c,alpha,beta,gamma, $ nfft1,nfft2,nfft3,order,ischrgd,verbose, $ checkacc,dsum_tol) implicit none integer nf ! _REAL_ a,b,c,alpha,beta,gamma,dsum_tol integer nfft1,nfft2,nfft3,order, $ ischrgd,verbose,checkacc c read unit cell parameters ! read(nf,*, end=10)a,b,c,alpha,beta,gamma c read the grid dimensions read(nf,*, end=11)nfft1,nfft2,nfft3,order, $ ischrgd,verbose,checkacc --- 3682,3698 ---- c------------------------------------------------------------------- c --- READ_EWALD --- ! subroutine read_ewald(nf,a,b,c,alpha,betta,gamma, $ nfft1,nfft2,nfft3,order,ischrgd,verbose, $ checkacc,dsum_tol) implicit none integer nf ! _REAL_ a,b,c,alpha,betta,gamma,dsum_tol integer nfft1,nfft2,nfft3,order, $ ischrgd,verbose,checkacc c read unit cell parameters ! read(nf,*, end=10)a,b,c,alpha,betta,gamma c read the grid dimensions read(nf,*, end=11)nfft1,nfft2,nfft3,order, $ ischrgd,verbose,checkacc *************** *** 3717,3723 **** 90 continue 100 continue ! call update_ucell_iso(a,b,c,alpha,beta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume, $ factor,verbose) return --- 3725,3731 ---- 90 continue 100 continue ! call update_ucell_iso(a,b,c,alpha,betta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume, $ factor,verbose) return *************** *** 4610,4631 **** c------------------------------------------------------------------- c --- UPDATE_UCELL_ISO --- c ...scales unit cell uniformly by factor; ! c alpha,beta,gamma unchanged... c ! subroutine update_ucell_iso(a,b,c,alpha,beta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume, $ factor,verbose) implicit none ! _REAL_ a,b,c,alpha,beta,gamma _REAL_ sphere,volume,factor(3) _REAL_ ucell(3,3),recip(3,3),dirlng(3),reclng(3) integer verbose _REAL_ facinv(3),result,distance integer i,j a = a * factor(1) b = b * factor(2) c = c * factor(3) volume = volume * factor(1)*factor(2)*factor(3) if ( verbose .eq. 1 ) write(6,99) a,b,c,volume 99 format(1x,'a,b,c,volume now equal to ',4f12.3) --- 4618,4643 ---- c------------------------------------------------------------------- c --- UPDATE_UCELL_ISO --- c ...scales unit cell uniformly by factor; ! c alpha,betta,gamma unchanged... c ! subroutine update_ucell_iso(a,b,c,alpha,betta,gamma, $ ucell,recip,dirlng,reclng,sphere,volume, $ factor,verbose) implicit none ! _REAL_ a,b,c,alpha,betta,gamma _REAL_ sphere,volume,factor(3) _REAL_ ucell(3,3),recip(3,3),dirlng(3),reclng(3) integer verbose _REAL_ facinv(3),result,distance integer i,j + # include "box.h" a = a * factor(1) b = b * factor(2) c = c * factor(3) + box(1) = a + box(2) = b + box(3) = c volume = volume * factor(1)*factor(2)*factor(3) if ( verbose .eq. 1 ) write(6,99) a,b,c,volume 99 format(1x,'a,b,c,volume now equal to ',4f12.3) ---------------------------------------------------------------------------- Make the following change to Makefile: ---------------------------------------------------------------------------- *** OLD Makefile --- NEW Makefile *************** *** 231,241 **** dynlib.o efield.o ene.o force.o indip.o mdread.o parallel.o period.o \ polder.o politr.o resnba.o runmd.o runmin.o sander.o sander_init.o \ ! setmm.o shake.o threeb.o tripl.o nonbon.o: box.h dpdynlib.o dpefield.o dpene.o dpforce.o dpindip.o dpmdread.o dpparallel.o \ dpperiod.o dppolder.o dppolitr.o dpresnba.o dprunmd.o dprunmin.o \ ! dpsander.o dpsetmm.o dpshake.o dpthreeb.o dptripl.o dpnonbon.o: box.h liblocmem.o libmdread.o librdparm.o librunmd.o librunmin.o libforce.o: box.h --- 231,242 ---- dynlib.o efield.o ene.o force.o indip.o mdread.o parallel.o period.o \ polder.o politr.o resnba.o runmd.o runmin.o sander.o sander_init.o \ ! setmm.o shake.o threeb.o tripl.o nonbon.o ewald.o: box.h dpdynlib.o dpefield.o dpene.o dpforce.o dpindip.o dpmdread.o dpparallel.o \ dpperiod.o dppolder.o dppolitr.o dpresnba.o dprunmd.o dprunmin.o \ ! dpsander.o dpsetmm.o dpshake.o dpthreeb.o dptripl.o dpnonbon.o \ ! dpewald.o: box.h liblocmem.o libmdread.o librdparm.o librunmd.o librunmin.o libforce.o: box.h ---------------------------------------------------------------------------- Temporary workarounds: none Routines affected: Ewald --