Boards:25 Issue

New Issue

polymerization setting
Shunsuke Mieda June 3, 2020, 9:14 a.m. #105

Dear all

Hello, I'm trying to run GRRM 17 program to examine polymerization reactions including Butadiene.

At first, I have tried to calculate the very simple reaction as below. However, any TS structure has NEVER been generated although GAMMA parameter was very large, SC-AFIR was changed to SC-AFIR2, or the different functionals were used.

To react these molecules, could anyone teach me what options should be tested, or how to set the input file.


S. Mieda

----Input file sample-----------


-1 1
C -6.04851 1.08683 -0.12952
C -4.65656 0.67303 -0.47198
C -3.56045 0.80445 0.37107
C -2.26684 0.40808 0.06896
H -6.39564 1.86874 -0.83809
H -4.49975 0.21350 -1.45756
H -3.73818 1.26439 1.35332
H -2.03118 -0.05733 -0.89259
C -7.06260 -0.07494 -0.17775
H -6.60275 -1.00092 0.30256
H -6.06390 1.54725 0.87543
H -1.44681 0.54598 0.77726
H -7.32889 -0.30211 -1.26268
C 3.12328 4.11991 -0.42148
C 1.91711 4.54010 0.00136
C 0.68735 4.38233 -0.75186
C -0.52022 4.79402 -0.32435
H 4.02317 4.26027 0.18117
H 3.24530 3.62228 -1.38871
H 1.82997 5.03655 0.97605
H 0.77488 3.88934 -1.72815
H -0.64143 5.28864 0.64465
H -1.42164 4.64671 -0.92286
H -8.00546 0.21979 0.39125
Add Interaction
GAMMA = 1000
Fragm.1 = 1-13,24
Fragm.2 = 14-23
1 2
GauMem = 2000
GauProc = 4

3 GRRM17 together with SIESTA
Eva M. Fernandez May 26, 2020, 6:43 a.m. #101

I'm trying to run GRRM17 together with SIESTA, so I've followed the instructions in the Manual web page (How to run GRRM with SIESTA).
I have correctly set all environment variables and created two sample files ( and h2o_min.inp).
But if I execute:
GRRM17p h2o_min
I get an error:
ioctl failed
What could be the cause?
Thanks in advance.
  • Re: GRRM17 together with SIESTA
    Kenichiro Saita

    In your environment (on your workstation), can you execute SIESTA as a standalone program (without GRRM17) without such errors?

    Please confirm your SIESTA program has been installed succesfully and you have set the environment variables "subsiesta" and "submpi" properly.

  • GRRM17 together with SIESTA
    Eva M. Fernandez

    Yes, the siesta program works Ok and I have set the environment variables "subsiesta" and "submpi" properly.

    Thanks in advance.

  • Re: GRRM17 together with SIESTA
    Kenichiro Saita

    Did you put the pseudopotential files (*.psf) in the same directory? In the sample case, two psf files O.psf and H.psf are required.

3 IRC Calculations
Cristina Nevado Feb. 15, 2020, 2:11 a.m. #96

Good afternoon,

I am working with GRRM17 to run IRC from transition states calculated in Gaussian.

When I use Gaussian to calculate IRC, I run 2 jobs, “reverse” and “forward”. Then, I can check that with “reverse” log file  I get the previous structure of the transition state and with “forward” file I can see the movement from the transition state to the product. When I run the IRC calculation in GRRM, I can see only the forward movement (from the transition state to the product).

Is it possible also to check the reverse movement from transition state to previous stratucture (or stating structure)?

Thank you very much

Best regards



  • Re: IRC Calculations
    Kenichiro Saita

    What options do you specify in your GRRM17 input file?

    In GRRM17, the IRC calculation always computes forward and backward IRC paths (unless the "Meta-IRC” option is specified).

    Of course, the forward and the backword paths can be computed from a saddle point (transition state). In other words, from a non-stationary point, only a mass-weighted steepest-descent path can be computed.

    Did you see the string "IRC FOLLOWING (FORWARD) STARTING FROM FIRST-ORDER SADDLE” in your GRRM17 log file? If you got "STEEPEST-DESCENT PATH FOLLOWING STARTING FROM NON-STATIONARY POINT”, it meant your input structure was not recognized as a saddle point.

    Even if your input structure was recognized as a 1st-order saddle point, it might correspond to the transition state structure between two conformers of the product.

    When you use the pre-optimized transition state structure by Gaussian (without GRRM17), such a problem sometimes appears. The SCF cycles might be converged in a different electronic structure. Try using the “MO Guess = filename.chk” option in order to read the Gaussian Checkpoint file as a proper initial guess.

    Or try the SADDLE calculation before the IRC calculation. The IRC paths can be automatically computed by specifying “Saddle+IRC” option.

    The DS-AFIR calculation would be able to help you to get the IRC paths between the reactant and the product. (Like the QST2 method in Gaussian.)

  • IRC Calculations
    Cristina Nevado

    Thank you very much for your kind reply.

    In my GRRM17 input file, I specify the option IRC. I see in my file.log "STEEPEST-DESCENT PATH FOLLOWING STARTING FROM NON-STATIONARY POINT" but at the end of the file also I see: "IRC following along both forward and backward directions were finished"

    I converted the file.log  in file.mol to see both movements in gaussview and, at this point, I can see only the forward movement. But as it is indicated in the file.log (IRC following along both forward and backward directions were finished), should doI see also the backward movement?

    Best regards


  • Re: IRC Calculations
    Kenichiro Saita

    No. The string "IRC following along both forward and backward directions were finished" just means the IRC calculation was normally terminated (termination message). Your GRRM17 output file (file.log) says that your input structure (pre-optimized TS structure obtained by Gaussian) did not converge in a saddle point in GRRM17. From a non-stationary point, only a mass-weighted steepest-descent path can be computed (the "backward" path cannot be theoretically defined).

    Try SADDLE calculation with "Saddle+IRC" and “MO Guess = filename.chk” options. You will get an actual TS structure and the "forward" and "backward" IRC paths.

    For example:

    # SADDLE/B3LYP/6-31G**

    0 1
    C  -0.079213112255    0.000038936045  -0.592031587948
    O   0.018800714554  -0.000018308130   0.717081381477
    H   0.322739833191   0.909186505065  -0.947006485140
    H   0.322739833191  -0.909359858691  -0.946723086445
    GauProc = 4
    GauMem = 100
    MO Guess = filename.chk

  • Re: IRC Calculations

    This issue seems to be solved by the above comment. (However, you can still add comments to this thread.)

1 Microiteration with external programs
Hiroko Satoh Oct. 18, 2019, 10:30 p.m. #93

Dear Developers,

I am interested in microiteration calculations with external programs (not with Gaussian). Could you please show us what data and format should be written in xxx_OUT4GEN.rrm. It seems that the manual has no descirptions about the case of TASK: MICROITERATION in xxx_INP4GEN.rrm.

I would appreciate if you could help us.

Best wishes,

Hiroko Satoh


  • The file xxx_OUT4GEN.rrm should be written in the same format in every case (TASK: "MAKE GUESS", "MICROITERATION", "ENERGY", "ENERGY and GRADIENT", or "ENERGY, GRADIENT, and HESSIAN").

    The required format, for example, for 3-atom system is:

    atom1  x1  y1  z1
    atom2  x2  y2  z2
    atom3  x3  y3  z3
    ENERGY =  E1  E2  E3
           =  E4  E5  E6
    S**2   =  <S2>
    DIPOLE =  μ(x)  μ(y)  μ(z)
      h(y1,x1)  h(y1,y1)
      h(z1,x1)  h(z1,y1)  h(z1,z1)
      h(x2,x1)  h(x2,y1)  h(x2,z1)  h(x2,x2)
      h(y2,x1)  h(y2,y1)  h(y2,z1)  h(y2,x2)  h(y2,y2)
      h(z2,x1)  h(z2,y1)  h(z2,z1)  h(z2,x2)  h(z2,y2)
      h(x3,x1)  h(x3,y1)  h(x3,z1)  h(x3,x2)  h(x3,y2)
      h(y3,x1)  h(y3,y1)  h(y3,z1)  h(y3,x2)  h(y3,y2)
      h(z3,x1)  h(z3,y1)  h(z3,z1)  h(z3,x2)  h(z3,y2)
      h(x3,z2)  h(x3,x3)
      h(y3,z2)  h(y3,x3)  h(y3,y3)
      h(z3,z2)  h(z3,x3)  h(z3,y3)  h(z3,z3)
      (x,x1)  (y,x1)  (z,x1)
      (x,y1)  (y,y1)  (z,y1)
      (x,z1)  (y,z1)  (z,z1)
      (x,x2)  (y,x2)  (z,x2)
      (x,y2)  (y,y2)  (z,y2)
      (x,z2)  (y,z2)  (z,z2)
      (x,x3)  (y,x3)  (z,x3)
      (x,y3)  (y,y3)  (z,y3)
      (x,z3)  (y,z3)  (z,z3)
      α(x,y)  α(y,y)
      α(x,z)  α(y,z)  α(z,z)

    If the file xxx_OUT4GEN.rrm is written in different formats, GRRM17 cannot read xxx_OUT4GEN.rrm. Therefore you need to put zero (0.0) into matrix elements of GRADIENT, DIPOLE, HESSIAN, DIPOLE DERIVATIVES, and/or POLARIZABILITY that cannot be obtained from your "sub link = aaa" program. Note: ENERGY requires 6 elements (E1E6). In usual cases, the energy from your "sub link = aaa" program have to be put into the first element (E1), and put zero into the others (E2E6). The other E2E6 elements are reserved for the multistate calculations or some special usages.

    In the case "TASK: MICROITERATION", GRRM17 requests you to perform microiterations by your "sub link = aaa" program and to give xxx_OUT4GEN.rrm which contains the coordinates after microiterations, energy, and gradients of the external atoms (although δV/δQn is almost zero for coordinates of external atoms Qn after microiterations). GRRM17 does not perform microiterations, and this is a specification (this might seem to be strange, but this is designed for general use).

    For further details, try a microiteration calculation with Gaussian program, and look at the xxx_LinkJOB.rrm file. You will understand how GRRM17 does a microiteration calculation. Even if you use the general interface, GRRM17 requires the same information.

    • Microiteration with external programs
      Hiroko Satoh

      Dear Prof. Saita,

      Thank you very much for the info and taking your time. We will give it a try.

      Best wishes,

      Hiroko Satoh


1 IRC Calculation
Cristina Nevado Sept. 16, 2019, 8:32 p.m. #83

I am interested to run IRC calculations with GRRM of some transition states that I obtained with Gaussian.

2 log files are generated in the calculation: name.log and name_GauJOB.log.

I run one file getting normal termination but when I open the file name_GauJOB.log in GaussView I cannot see the expected movement of the molecule (reverse or forward).

Maybe I did not indicate in the file the correct parameters. 

My file:

 # IRC/b3lyp/GenECP

0 1


GauMem = 200
GauProc = 8
Bond Condition
12 75 < 2.06
Add Interaction
1 2

GAMMA = 200

C H N O F Si

Au 0

 In this example, I indicate as Bond Condition 12 75: I am working with transition states:

 Here, do I have to indicate the numbers of the 2 atoms involved in the formation of the new bond, right?

  1. When I open the log file of my transition state (obtanied in Gaussian) in gaussview I see the movement of the atoms when I select the option Vibrations. Which is the distance that I have to indicate for the calculation in GRRM?. Is the distance that I have before to select vibrations? The distance that I have to indicate is the difference between the distance of these atoms in the product and the precursor? or between the transition state and the product?
  2. Do I indicate > or < based on what?

 Thank you very much in advance for your kind help


  • Re: IRC Calculation
    Kenichiro Saita

    The GRRM17 program uses Gaussian program to only compute single point calculation, and the file name_GauJOB.log is the Gaussian log file of the last step of your calculation. In this case, it would correspond to the freq calculation log file at the minimum point of the backward IRC path, so that you cannot see the IRC path information in the file name_GauJOB.log with GaussView. The IRC information is stored in the file name.log. To visualize the structure in IRC path computed by GRRM17 with GaussView, you need to plot the corresponding structure and energy separately from the file name.log.

    In the "Bond Condition ... END" option, "12 75 < 2.06" means "the distance between atoms 12 and 75 to be less than or equal to 2.06 angstroms". This bond condition will be appled to the SC-AFIR search, but not to the IRC path calculation. Do not confuse this option with the IRC=Phase=(N1,N2 [,N3 [,N4]]) option in Gaussian. In GRRM17, an IRC calculation always computes forward and backward IRC paths.

    The options "NoBondRearrange", "Bond Condition ... END", and "Add Interaction ... END" have no effect on the IRC calculation, because "IRC" keyword just performs the IRC path calculation. Such options are effective in the "SC-AFIR" calculation. In SC-AFIR, the "Bond Condition ... END" option cannot be used simultaniously with the "NoBondRearrange" option (the latter keyword in the input file is only effective).

    You may need to add the basis to the valence part of Au:

    C H N O F Si 0
    Au 0

    Au 0

  • Re: IRC Calculation

    This issue seems to be solved by the above comment. (However, you can still add comments to this thread.)

AFIR Web unavailable: 31 Aug - 2 Sep 2019
afiradmin Aug. 20, 2019, 1:41 p.m. #79

Due to a planned power outage in Hokkaido University, all services on our website will be unavailable during the following date and times:

31 August 2019, 15:00 (GMT+9) - 2 September 2019, 13:00 (GMT+9)

Thank you for your cooperation.

  • AFIR Web server has been back

    All services on our website are running now.

    Thank you for your cooperation.

Elementary questions
YU Kaneko June 28, 2019, 9:55 a.m. #76

Dear GRRM developers,

Very recently, I became a new user of GRRM17. I managed to perform MC-AFIR calculation for exploring new reaction between different two molecules.  I have three question now.

Question 1 :  I would like to understand the meaning of "Fragment" and "Part". In "MC-AFIR", "part" means the index of molecules for generating initial configration of geometry. I found same question and answer (  Applying artificial force in SC-AFIR , "Fragment" means the index of atoms that is applied to AFIR external force during reaction searching.   If so, if I could guess the atoms of reaction area properly, the whole calculation time could be reduced.

Question 2: What is the usage of "GRRM17.out" ?  If it is something for summarizing reaction pathway, I would like to use it.

Question 3: When I could not found the desired or supposed reaction pathway among the extracted reaction pathways, could  large value setting for NFault or GAMMA be effective ? In my understanding, NFault just limits the number of reaction patth collection, and GAMMA limits the walkable energy hieght during TS survey. If you have effective option,  I would like to know.

Best regards,

Yu Kaneko 

GRRM linked with MOLPRO2015
Yasuhiro Shigemitsu June 25, 2019, 11:45 a.m. #75

Dear GRRM developers,

I'm trying to use GRRM17 combined with MOLPRO2015 but the combination does not work.  The GRRM tutorial reads that the version 2006, 2008,2010,2012 were successfully tested with GRRM. Kind helps or information would appreciated to use the combination, GRRM17 plus MOLPRO2015.

Best regards,

Yasuhiro Shigemitsu

  • GRRM17 linked with Molpro 2015 or later
    Kenichiro Saita

    GRRM17 does not recognize the keyword "%link=molpro2015" or later.  Consequently, GRRM17 is officially compatible with Molpro 2006, 2008, 2010, and 2012. However, you can use Molpro 2015 or later at your own risk.


    With Molpro 2015: 

    Whereas we have not fully tested Molpro2015 jobs, the formats of output files of Molpro 2015 seems to be the same as those of Molpro2012, and GRRM17 will work with Molpro2015 by using the link keyword "%link=molpro2012". 


    With Molpro 2018 or later: 

    GRRM17 knows a Molpro job has been terminated normally if the string "Variable memory released" is appeared in the Molpro output file (XXX_MolJOB.out), because Molpro2015 or earlier prints the string in the last line of the output. But, Molpro2018 or later uses "Molpro calculation terminated" instead of "Variable memory released". Thus, GRRM17 will never understand the termination of the Molpro2018 (or later) calculation, and then GRRM17 program will stop with an error.  This issue will be evaded by using the TEXT card in Molpro input. For example: 

    Molpro Input



    $STR='Variable memory'
    TEXT,$STR released
    RSPT2 STATE 2.1

    We are aware of another issue; the format of the punch file (XXX_moljob.pun) is different from that of Molpro2015 or earlier, then GRRM17 cannot read the energies from the punch file. This issue will be evaded by setting "submol" command to the shellscript in which the punch file format to be converted. For example: 

       setenv submol "./script.csh"   (csh/tcsh case)

    and prepare "script.csh" like as 

    #!/bin/csh -f
    set PUN=`echo $2:r.pun | tr \[A-Z\] \[a-z\]`
    set TMP='tmp'

    molpro -W./ -d./ $1 $2

    grep "Energy" ${PUN} > ${TMP}
    set i = 1
    set n = `wc -l ${TMP} | awk '{print $1}'`
    while ( $i <= $n )
       set line = `head -n $i ${TMP} | tail -n 1`
       echo $line >> ${PUN}
       @ i = $i + 1
    \rm -f ${TMP}
    unset i, n, line
    unset PUN, TMP


    Alternatively, you can use the general interface in GRRM17. Please see "General interface with external ab initio programs".

    We are aware of another issue; the format of the punch file (XXX_moljob.pun) is bit different from that of Molpro2015 or earlier, then GRRM17 cannot read the energies from the punch file. For example:

    RSPT2 STATE 1.1 Energy    -114.21130105   (Molpro2015 or earlier)

    RSPT2 STATE  1.1 Energy    -114.21130105   (Molpro2018 or later)

    So, make sure whether you provide an appropriate string to the Molpro Input section. For example,

    Molpro Input

    $STR='Variable memory'
    TEXT,$STR released
    RSPT2 STATE  2.1


    We have made several tests using Molpro 2015.1 (patch level 33), Molpro 2018.2, and Molpro 2019.1.2, but we cannot give any warranty.

Taiki Kato Nov. 27, 2018, 10:40 a.m. #71

Dear Prof. Maeda & laboratory people,
I'm trying ONIOM method on MC-AFIR.
But, I don't understand the exact format of the ONIOM layer & MC-AFIR part.
This is my input file with error calculation log.


# MC-AFIR/ONIOM(uwb97xd/6-31G : upm6) integral(coarsegrid)
0 1
H  -5.09198232  1.46148646  0.64415956 H   1
H  -5.05103711  1.51846702 -1.80580446 H   1
H  -5.08619301  3.61195144 -0.53172902 H   1
H  -3.07552771  2.21038282 -0.53072587 H   1
Si -4.57617575  2.20057460 -0.55601978 H   1
Si  0.92262713 -0.42951851 -0.06586408 H   2
H   0.32586471  0.89377867  0.44846657 H   2
H   0.92141623 -0.42839950 -1.60585939 H   2
Si -0.36596427 -2.22060001  0.71381911 H   2
Si  0.54148635 -4.23086015 -0.06762047 H   2


Si  7.89452207  2.71361008 -0.06626650 L H 16 2
Si  6.60446476  0.92311796  3.05275285 L H 16 2
Si  3.11867252 -0.64788568  3.05254379 L H 12 2
Si -2.56235026 -2.00324774 -0.06393534 L H  9 2
Si  1.83058292 -2.43960890  3.83181251 L H 19 2
Si  1.45021923 -6.24149307  3.83086153 L H 19 2
Si -0.74565359 -6.02268238  0.71234080 L H 10 2
Si  2.35719613 -8.25146150  3.04887994 L H 23 2
Si  4.55185684 -8.47158271  3.82932265 L H 33 2
External Atoms
Si  7.50985596 -1.08833081  6.95223274 L   2
Si  7.12790465 -4.89067932  6.95082004 L   2
Si  8.03487102 -6.90152356  6.17032595 L   2
Si  8.41640895 -3.09947186  6.17176893 L   2


H   6.74535629 -8.69308650  8.48930406 L   2
H   7.50965005 -1.08890581  8.49222140 L   2
H   4.02610987 -2.65912846  8.49169783 L   2
H   3.64455597 -6.46260745  8.48999760 L   2
H   0.54408428 -4.23167248  8.49284240 L   2
MicroIt = (MMOnly)
GauProc = 4
NSample = 20
MinFC = -1
MC = ReactivePathOnly
Add Interaction
Fragm.1 = 1-5
Fragm.2 = 6-43
Fragm.3 = 66-145
1 2
1 3 -
GAMMA = 600


I think that the number order of ONIOM layer or MC-AFIR part is not good.
Please teach me the exact number order of ONIOM layer & MC-AFIR part.
Best regards

  • Re: ONIOM on MC-AFIR
    Satoshi Maeda

    To combine ONIOM and MC-AFIR, please use the following format.

      atom  x  y  z  part-number  layer  [link-atom info]


      atom :  atom type
      x :  X-coordinate
      y :  Y-coordinate
      z :  Z-coordinate
      part-number :  part number for MC-AFIR (see Part designation for random structure generation)
      layer :  layer assignment (H/M/L) for ONIOM
      [link-atom info] :  (as needed) link-atom parameters for ONIOM

    In your case, for example:

    H  -5.09198232  1.46148646  0.64415956  1  H
    H  -5.05103711  1.51846702 -1.80580446  1  H
    H  -5.08619301  3.61195144 -0.53172902  1  H
    H  -3.07552771  2.21038282 -0.53072587  1  H
    Si -4.57617575  2.20057460 -0.55601978  1  H
    Si  0.92262713 -0.42951851 -0.06586408  2  H
    H   0.32586471  0.89377867  0.44846657  2  H
    H   0.92141623 -0.42839950 -1.60585939  2  H
    Si -0.36596427 -2.22060001  0.71381911  2  H
    Si  0.54148635 -4.23086015 -0.06762047  2  H
    Si  7.89452207  2.71361008 -0.06626650  2  L H 16
    Si  6.60446476  0.92311796  3.05275285  2  L H 16
    Si  3.11867252 -0.64788568  3.05254379  2  L H 12
    Si -2.56235026 -2.00324774 -0.06393534  2  L H  9
    Si  1.83058292 -2.43960890  3.83181251  2  L H 19
    Si  1.45021923 -6.24149307  3.83086153  2  L H 19
    Si -0.74565359 -6.02268238  0.71234080  2  L H 10
    Si  2.35719613 -8.25146150  3.04887994  2  L H 23
    Si  4.55185684 -8.47158271  3.82932265  2  L H 33

2 Metal Element Containing Surface Reaction Analysis
Taiki Kato Nov. 8, 2018, 5:30 p.m. #67

Dear Prof. Maeda & Laboratory People
I'm struggling to analyze metal element containing surface reaction analysis such as { TiCl4 on TiN }, { EtOH on HfO2 }, { SiH4 on Ru } and so on.
I use GRRM17 which is not feasible PBC surface model.
Then, I use cluster surface model on these metal element containing surfaces.
However, it is difficult for me to converge these surfaces { TiN, HfO2, Ru } SCF calculation by Gaussian.
I also use Materials Studio DMol3 which is easy to converge SCF calculation but not suitable for GRRM.
Please teach me how to converge these SCF calculation by Gaussian or some alternative convergent technology.

  • Re:
    Satoshi Maeda

    Although GRRM17 doesn't contain an interface with DMol3, GRRM17 can be used together with any energy computation code if a user prepare a simple interface script. Please see the following page:

    • Other QM code trial
      Taiki Kato

      Thank you for your kind reply.

      I'll consider the GRRM application trial not only on Gaussian but also on the other QM codes.

      DMol3 doesn't have any analytical hessian calculation method.

      Therefore, I think that DMol3 may be very slow on GRRM.

      And I hope to choice more established way because I'm a biginner on these operations.

      How is the performance of SIESTA (or Turbomole) on cluster model analysis ?

      Will you please teach me SIESTA (or Turbomole) setting manuals ?

      If you know good QM-codes, please give me some advice.

    • Re: Other QM code trial
      Kenichiro Saita

      You can try SIESTA program. Please see the manual page "How to run GRRM with SIESTA."

      NOTE: Though SIESTA 3.2 and older versions are distributed free of charge to only academics, SIESTA 4.0 and later vaesions are released under the terms of the GPL open-source license.

      Turbomole program is also available in GRRM17 (see How to run GRRM with TURBOMOLE) in terms of cluster model analyses. The RI-approximation would be effective in computational time.

      DFT calculations for periodic systems are available in Turbomole (7.2 and later), but current version of GRRM17 binaries does not support the RIPER module of Turbomole. If you use Turbomole (7.2 or later) and its RIPER module, use General interface at your own responsibility.