I have tested a SC-AFIR calculation using your sample input; the artificial force was not applied between any fragments because the two moieties were too far away (the input structure (EQ0) had judged as a dissociation channel (DC) structure by the default parameter DownDC=8, and you would be able to see the string "Not to be applied" in the file jobname_EQ_test.rrm.). Therefore, no reaction paths had been computed and no TS structures had also been generated.

Does the intermolecular distance have a meaning (e.g. experimentally solved structure)? If no, it might be worth to try a MC-AFIR calculation. In addition, I think "Stable=Opt" option is recommended in your case.


----- Sample Input for MC-AFIR -----

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

--------------------

This is just my comment, and this might not be a solution to your problem...

Saita-Sensei

Thank you for the great advice!

I'm trying to shorten the distance between two fragments, to use the options, and to generate multi-structires by using MC-AFIR.

BTW... could you teach me why the fragment ranges were changed? Did the wide ranges have bad influence on the reaction?

Best regards,

S. Mieda

In the above sample input for MC-AFIR, only the sp² carbon atoms are selected for the fragment. It would be a natural choice in terms of the chemical knowledge, but it was not a strong suggestion. In order to make a global search, the wider fragment range would be required, but it would take a long computational time. So, I chose only the sp² carbon atoms to the fragments as a sample job.

NOTE: The terms "part" and "fragment" are used in different meanings in GRRM17 program:

· "Fragment" means the fragment to which the artificial force is applied by the AFIR function.
· "Part" means a moiety of the system (or a monomer in the molecular cluster). The defined parts to be used for random geometry generation (all parts to be randomly distributed).

The fragment is not necessarily identical to the part. For details, please see the manual page.
In the SC-AFIR calculation, the fragments to be automatically defined, so that you do not have to put the lines "Fragm.1 = 1-13,24" and "Fragm.2 = 14-23" into the SC-AFIR input file.

(For your information) this might be an example input file for SC-AFIR local search:

----- Sample Input for SC-AFIR -----

# SC-AFIR/B3LYP/6-31G
 
-1 1
C -6.04851  1.08683 -0.12952  1
C -4.65656  0.67303 -0.47198  1
C -3.56045  0.80445  0.37107  1
C -2.26684  0.40808  0.06896  1
H -6.39564  1.86874 -0.83809  1
H -4.49975  0.21350 -1.45756  1
H -3.73818  1.26439  1.35332  1
H -2.03118 -0.05733 -0.89259  1
C -7.06260 -0.07494 -0.17775  1
H -6.60275 -1.00092  0.30256  1
H -6.06390  1.54725  0.87543  1
H -1.44681  0.54598  0.77726  1
H -7.32889 -0.30211 -1.26268  1
C  3.12328  4.11991 -0.42148  2
C  1.91711  4.54010  0.00136  2
C  0.68735  4.38233 -0.75186  2
C -0.52022  4.79402 -0.32435  2
H  4.02317  4.26027  0.18117  2
H  3.24530  3.62228 -1.38871  2
H  1.82997  5.03655  0.97605  2
H  0.77488  3.88934 -1.72815  2
H -0.64143  5.28864  0.64465  2
H -1.42164  4.64671 -0.92286  2
H -8.00546  0.21979  0.39125  1
Options
NRUN = 1
Add Interaction
Target = 2-4,14-17
GAMMA = 200
END
SC = InterOnly
Stable = Opt
GauMem = 2000
GauProc = 4
 

--------------------

NRUN is used in order to generate a random initial structure.
SC = InterOnly limits the number of paths to be computed (artificial force to be applied between different parts only).

This is just an example. The SC-AFIR has various options, so you may wish to experiment to decide which one fits your needs best.

Saita-sensei

Thank you very much for the useful comments.

I understand the difference between Frangments and Parts (and Targets).

Especially, SC = InterOnly option seems to be very useful for my calculation. I'm trying to use the other options for SC-AFIR, too.

If the reaction is carried out, I write the input file here to share the information.

Best regards,

S. Mieda 

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.

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

Thanks in advance.

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.

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.)

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

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
Options
Saddle+IRC
GauProc = 4
GauMem = 100
MO Guess = filename.chk

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

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:

RESULTS
CURRENT COORDINATE
atom1  x₁  y₁  z₁
atom2  x₂  y₂  z₂
atom3  x₃  y₃  z₃
ENERGY =  E₁  E₂  E₃
       =  E₄  E₅  E₆
S**2   =  <S²>
GRADIENT
  g(x₁)
  g(y₁)
  g(z₁)
  g(x₂)
  g(y₂)
  g(z₂)
  g(x₃)
  g(y₃)
  g(z₃)
DIPOLE =  μ(x)  μ(y)  μ(z)
HESSIAN
  h(x₁,x₁)
  h(y₁,x₁)  h(y₁,y₁)
  h(z₁,x₁)  h(z₁,y₁)  h(z₁,z₁)
  h(x₂,x₁)  h(x₂,y₁)  h(x₂,z₁)  h(x₂,x₂)
  h(y₂,x₁)  h(y₂,y₁)  h(y₂,z₁)  h(y₂,x₂)  h(y₂,y₂)
  h(z₂,x₁)  h(z₂,y₁)  h(z₂,z₁)  h(z₂,x₂)  h(z₂,y₂)
  h(x₃,x₁)  h(x₃,y₁)  h(x₃,z₁)  h(x₃,x₂)  h(x₃,y₂)
  h(y₃,x₁)  h(y₃,y₁)  h(y₃,z₁)  h(y₃,x₂)  h(y₃,y₂)
  h(z₃,x₁)  h(z₃,y₁)  h(z₃,z₁)  h(z₃,x₂)  h(z₃,y₂)
  h(z₂,z₂)
  h(x₃,z₂)  h(x₃,x₃)
  h(y₃,z₂)  h(y₃,x₃)  h(y₃,y₃)
  h(z₃,z₂)  h(z₃,x₃)  h(z₃,y₃)  h(z₃,z₃)
DIPOLE DERIVATIVES
  dμ(x,x₁)  dμ(y,x₁)  dμ(z,x₁)
  dμ(x,y₁)  dμ(y,y₁)  dμ(z,y₁)
  dμ(x,z₁)  dμ(y,z₁)  dμ(z,z₁)
  dμ(x,x₂)  dμ(y,x₂)  dμ(z,x₂)
  dμ(x,y₂)  dμ(y,y₂)  dμ(z,y₂)
  dμ(x,z₂)  dμ(y,z₂)  dμ(z,z₂)
  dμ(x,x₃)  dμ(y,x₃)  dμ(z,x₃)
  dμ(x,y₃)  dμ(y,y₃)  dμ(z,y₃)
  dμ(x,z₃)  dμ(y,z₃)  dμ(z,z₃)
POLARIZABILITY
  α(x,x)
  α(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 (E₁ – E₆). In usual cases, the energy from your "sub link = aaa" program have to be put into the first element (E₁), and put zero into the others (E₂ – E₆). The other E₂ – E₆ 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/δQ_n is almost zero for coordinates of external atoms Q_n 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.

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:

GauInpB
C H N O F Si 0
6-31g**
****
Au 0
SDD
****

Au 0
SDD
END

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

All services on our website are running now.

Thank you for your cooperation.

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
BASIS=6-311+G(d,p)

{MCSCF,maxit=39,GRADIENT=1.d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,12;
Closed,4;
WF,16,1,0;
STATE,2;}

{RS2,Mix=2,Root=2,SHIFT=0.3,MaxIt=99,MaxIti=999;
ORBITAL,IGNORE_ERROR;
WF,16,1,0;
STATE,2;}

FORCE;
IF(STATUS.LT.0) STOP
$STR='Variable memory'
TEXT,$STR released
---
RSPT2 STATE 2.1
END

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
end
\rm -f ${TMP}
unset i, n, line
unset PUN, TMP

exit

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
︙
IF(STATUS.LT.0) STOP
$STR='Variable memory'
TEXT,$STR released
---
RSPT2 STATE  2.1
END

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.

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

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

where

  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
...

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:

https://afir.sci.hokudai.ac.jp/documents/manual/48