qub (c46b2)

Syntax for QUB Command

QUB { QCP } [atom-selection] [other-spec] [unit-spec] [qrxn-spec]
{ BQCP bqcp-spec }
{ SQCP }
{ BFEP bqcp-spec }
{ SFEP }
{ QCOP }

bqcp-spec::= [KLEVel int]
other-spec::= [BEADs int] [NBMOve int] [MCONfiguration int]
[MCEQquilibration int] [TEMP real] [QRXN qrxn-spec]
[TIAC] [CHAC] [NOEW] [IRAN] [FAST] [FFOCK]
unit-spec::= [FIRST int] [OUNI int] [OUNJ int] [BDIN int] [BDOUT int]
[NUNI int] [BEGI int] [STOP int] [SKIP int]
qrxn-spec::= [RXNA int] [RXNB int] [RXNC int] [DELD real] [RESL real]

Keyword Deflt Purpose
------- ----- ----------------------------------------------------------
QCP off Use Quantized Classical Path method with Metropolis sampling.
BQCP off Use Quantized Classical Path method with bisection sampling.
KLEV must be specified.
SQCP off Use Quantized Classical Path method with staging sampling.
BFEP off Use PI-FEP/UM which entails BQCP with mass-perturbation.
SFEP off Use PI-FEP/UM which entails SQCP with mass-perturbation.
QCOP off Use open path-integrals with staging algorithm.
BEADS 32 Number of quasiparticles (beads) per classical particle.
NBMO 1 Number of beads to be changed in each MC move. Must be 1
when using Bisection algorithm. When using Metropolis sampling
must be greater than 1 (due to centroid approximation).
KLEV 5 Should be specified when using Bisection algorithm. In general,
BEADS => 2**KLEV. Optimal sampling is achieved when
BEADS = 2**KLEV,
i.e. BEADS 8 KLEV 3, BEADS 16 KLEV 4, BEADS 32 KLEV 5, etc.
TIAC off Use Takahashi-Imada higher-order correction
CHAC off Use Chin higher-order correction
FAST F Use fast energy routines (avoids gradient calculation). Only
QUANTUM module.
FFOCK F Use fast Fock matrix updating. Only QUANTUM module.
NOEW F Do not use Ewald summation
MCON 10 Number of MC moves for averaging
MCEQ 10 Number of MC moves for equilibration
TEMP 298.15 Temperature (should be same as used for classical trajectory)
QRXN off Use reaction coordinate analysis
RXNA 99999 PSF number of atom A
RXNB 99999 PSF number of atom B
RXNC 99999 PSF number of atom C; if excluded reaction coordinate is A-B
DELD 0.01 Bin width along reaction coordinate
RESL 10.0 Resolution in bin analysis. Rarely needs to be adjusted.
FIRS -1 Unit number of first trajectory file
OUNI 6 File unit for QUB output (Qcl->qm. Partition function qm
correction)
OUNJ 6 File unit for QUB output for second isotope (w/QFEP)
OUNK 6 File unit for end-to-end distribution in open path-integral simulation
OUNL 6 File unit for end-to-end distribution, second isotope
BDIN 0 File number for coordinates of all beads (from a previous job)
BDOU 0 File number for coordinates of all beads to be written out
at the end of the current job
NUNI 1 Numbero of trajectory files, beginning with FIRS
BEGI 0 Beginning step number of trajectory file
(see READ TRAJECTORY FOR DETAIL)
STOP 0 See read traj
SKIP 1 See read traj
IRAN -1 Random number generator seed, default value switches on time
dependent random seed (recommended). For runs on the same
trajectory, but with different isotopes, the same random
number should be used (e.g. use -1 for "isotope 1" run and
the random number from the output of that run (printed at
beginning of QUB output) for "isotope 2" run.


Top
Assuming the following CHARMM variables have been set:
! Path-integral settings
set NSkip 1 !
set NSave 100 !
calc NSkip @NSave*@NSkip ! Skip # of configurations from trj
set NBeads 16
set T 298.15
set Neq 100 ! Equilibration of beads not necessary
! when bisecting entire bead
set Nav 100 ! Use ~ 1000
set NMove 1
set KLev 4 ! 2**KLev should equal number of beads
set RNum -1 ! Time-dependent randum number seed

1) Add average QM effect

...

! read coordinates from trj file
open read unit 11 unfo name @scrDIR/@MOD.@d0.@j.trj
open write unit 12 form name @resDIR/@MOD.@d0.N@NBeads.qub ! QUB output

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Run PI
QUB BQCP SELE QMS .AND. (TYPE C5 .OR. TYPE H7 .OR. TYPE O16) SHOW END
FAST FFOCK -
KLEV @KLev TEMP @T MCON @Nav MCEQ @Neq BEAD @NBeads NBMOVE @NMove -
IRAN @RNum -
FIRST 11 OUNI 12 BDIN -13 BDOUT -14 -
NUNI 1 BEGI -1 STOP -1 SKIP @NSkip

...


2) Add QM effect along a reaction coordinate of the type A-B---C, where B is
transferred from A to C. Typically, this would be done seperately for each
window of an umbrella sampling.
If atom C is not defined in the input, the reaction coordinate is taken as
the A-B distance.
Typically, several different runs are done along a reaction coordinate and
the results (from *.qub output files) combined using a post-processing
program (available upon request) to obtain average along reaction
coordinate.

...

! read coordinates from trj file
open read unit 11 unfo name @scrDIR/@MOD.@d0.@j.trj
open write unit 12 form name @resDIR/@MOD.@d0.N@NBeads.qub ! QUB output

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Run PI
! Reaction coordinate is defined as A-B---C
! C5(5)-H7(7)---O16(16)
QUB BQCP SELE QMS .AND. (TYPE C5 .OR. TYPE H7 .OR. TYPE O16) SHOW END
FAST FFOCK -
KLEV @KLev TEMP @T MCON @Nav MCEQ @Neq BEAD @NBeads NBMOVE @NMove -
QRXN RXNA 5 RXNB 7 RXNC 16 DELD 0.01 RESL 10.0 -
IRAN @RNum -
FIRST 11 OUNI 12 BDIN -13 BDOUT -14 -
NUNI 1 BEGI -1 STOP -1 SKIP @NSkip

...

3) Same as 2, but using different isotope.
Add the following lines prior to calling
the QUB command. Everything else is identical to exmape 2).

...

! Change mass of relevant atom
set mass 2.0140

scalar mass set @mass select QMS .AND. TYPE H7 show end

! Section 2) comes here

4) Mass perturbation.

...

! read coordinates from trj file
open read unit 11 unfo name @scrDIR/@MOD.@d0.@j.trj
open write unit 12 form name @resDIR/@MOD.@d0.N@NBeads.h.qub ! isotope 1 QUB output
open write unit 13 form name @resDIR/@MOD.@d0.N@NBeads.d.qub ! isotope 2 QUB output

! D isotope
! Change mass of relevant atom in WMAIN array
set dmass 2.0140 ! Assign mass variable
scalar 1 = wmain
scalar wmain show
scalar wmain = mass ! Assign mass array to wmain
scalar wmain show
scalar wmain set @dmass select QMS .AND. TYPE H7 show end
scalar wmain show

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Run PI
! Reaction coordinate is defined as A-B---C
! C5(5)-H7(7)---O16(16)
QUB BFEP SELE QMS .AND. (TYPE C5 .OR. TYPE H7 .OR. TYPE O16) SHOW END
FAST FFOCK -
KLEV @KLev TEMP @T MCON @Nav MCEQ @Neq BEAD @NBeads NBMOVE @NMove -
QRXN RXNA 5 RXNB 7 RXNC 16 DELD 0.01 RESL 10.0 -
IRAN @RNum -
FIRST 11 OUNI 12 OUNJ 13 BDIN -14 BDOUT -15 -
NUNI 1 BEGI -1 STOP -1 SKIP @NSkip

! Reset wmain
scalar wmain = 1
scalar wmain show

...