# qub (c49b1)

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.

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

...

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

...