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
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
COORDINATES
Options
GauMem = 200
GauProc = 8
NoBondRearrange
Bond Condition
12 75 < 2.06
Add Interaction
Fragm.1=1-40
Fragm.2=41-75
1 2
GAMMA = 200
END
GauInpB
C H N O F Si
6-31g**
****
Au 0
SDD
END
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?
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:
GauInpB
C H N O F Si 0
6-31g**
****
Au 0
SDD
****
Au 0
SDD
END
Re: IRC Calculation
afiradmin
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
afiradmin
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 , https://afir.sci.hokudai.ac.jp/documents/manual/32) "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
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.
ONIOM on MC-AFIR
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
Options
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
END
-------------------------------------------------------------------------------------
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]
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
...
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.
Sencerely
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.
1
SC-AFIR option 'SC=NegativeON'
Ayako Furuhama
Aug. 8, 2018, 7:25 a.m.
#64
Dear Professor Maeda and his group members;
SC-AFIR calculation is performing to obtain various initial reactants of a reaction. But the following calculation is not stopped. [I do not know the reasons... Maybe 'SC=NegativeON' (?) This was used in previous calculations.] In this case, I will not follow PT structures and their TS structures due to computational cost.
Please let us know the additional or unnecessary options for such calculations. input_PARAM.rrm is following:
INPUT DATA SET OF THE GRRM VER. 13.x
-------------------------------------------------------------------------
Energy Calculation = GAUSSIAN03
-------------------------------------------------------------------------
-------------------------------------------------------------------------
Job type = 16
-------------------------------------------------------------------------
UB3LYP/6-31G*
-------------------------------------------------------------------------
0 1 XX XX 0
-------------------------------------------------------------------------
XYZ coordinate
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-------------------------------------------------------------------------
NRun = 32
IniStr = Random
GauMem = XXX (MW)
GauProc = X
ADD INTERACTION
SC = NegativeON
DownDC = 15
RandDC = 15
STABLE = OPT
READ BOND CONDITION
Thank you very much!
Option for SC-AFIR
Ayako Furuhama
Dear Prof. Maeda,
I am sorry but I clarify my question.
I search TS and EQ structures for the initial reaction of A + B system by using SC-AFIR/UB3LYP/6-31G(p).
A is a tri-atomic molecule and B is an organic molecule including one nitrogen and some carbon and hydrogen atoms.
After one month starting a SC-AFIR calculation (B= one nitrogen + three carbons with some hydorgens) with options (i.e. NRUN=32,GauMem=1000,GauProc=4,Add Interaction,and Gamma=100 alog with SC=NegativeON,SC=LUP-ON,NoBondRearrange,DownDC=15,RandDC=15,Stable=OPT), 343EQ and 629TS (as well as 5900 TS in PT_list) were found. But this job does not done yet. EQ and TS searches are not stopping.
Could you please inform how to finish the SC-AFIR calculation properly? Should I perporm FirstOnly calculations in advance to found out proper initial conformation?
Thank you very much in advance.
Ayako Furuhama
Re:
Satoshi Maeda
I don't know your system well, but I guess that many similar (but slightly different) EQs are generated in the flat regions of the PES in your calculation. This often happens when a spin unrestricted DFT method is used as an energy calculation method.
One can avoid this situation by one of the following three ways.
1. By increasing the grid density in DFT calculations.
E.g.) UB3LYP/6-31G(d) Int(Grid=99590)
This would be the most straightforward way, but DFT calculations with dense grids will be computational demanding.
2. By using a structural clustering option.
E.g.) MatchDecScale = 3.0
See the option list about this option: https://afir.sci.hokudai.ac.jp/documents/manual/51
3. By applying a weak force between all atom pairs.
E.g.) A simple case for CO+OH is shown below:
# SC-AFIR/UB3LYP/D95V
0 2
C 0.0 0.0 0.0 1
O 1.2 0.0 0.0 1
O 3.0 0.0 0.0 2
H 4.0 0.0 0.0 2
Options
Add Interaction
Fragm.1=1
Fragm.2=2
Fragm.3=3
Fragm.4=4
1 2 3.33
1 3 3.33
1 4 3.33
2 3 3.33
2 4 3.33
3 4 3.33
Gamma=200.0
END
EQOnly
KeepLUPPATH
With the above input, a weak force with Gamma=3.33 kJ/mol is added to all atom pairs during the SC-AFIR calculation with Gamma=200.0 kJ/mol. A recommended value of the weak force is 10.0/(N-1) kJ/mol, where N is the number of atoms in the system. All EQs and TSs obtained are those on the AFIR function consisting of these force terms and thus are not actual ones. Therefore, after the initial SC-AFIR calculation, one needs to perform a RePATH calculation reading the results of the initial SC-AFIR calculation to eliminate the weak force. Since TSs obtained by the initial SC-AFIR calculation are not actual TSs, it is recommended to avoid computing TSs accurately during the initial search using the EQOnly and KeepLUPPATH options (also see ReadBareEnergy option).
When the weak force is added, local minima for weak complexes between two molecules are formed in regions where the two molecules have sufficient interactions. When DFT grids that are not dense enough are used in spin unrestricted DFT calculations, many imaginably local minima can be created in the flat regions of the PES where two molecules are arranged far away from each other. The weak force help avoids this problem. I emphasize again that the weak force must be eliminated afterward by the additional RePATH calculation.
An example input for the RePATH calculation is also shown below:
%infile=xxx
# RePATH/UB3LYP/D95V
0 2
C 0.0 0.0 0.0 1
O 1.2 0.0 0.0 1
O 3.0 0.0 0.0 2
H 4.0 0.0 0.0 2
DS-AFIR issue 2
Seiji Mori
July 27, 2018, 11:33 a.m.
#57
Using the DS-AFIR to calculate the path that connects two structures, we found out that the calculation showed a final message as follows:
The structure is dissociating
OPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTOPTO
Path-MIN-optimization failed.....
If the error is caused because one of the intermediates is dissociating, we are wondering how we can solve this issue?
I tried using the SC=FindDC along with %infile=NAME command to retrieve the information obtained in previous calculations but I obtained the same error. Is there a way to avoid the failure due to structure dissociation so the DS-AFIR is completed?
Thank you very much in advance.
Re:
Satoshi Maeda
DS-AFIR issue
Seiji Mori
July 12, 2018, 3:19 p.m.
#55
During a DS-AFIR calculation, AFIR program found 5 Equilibrium Structures (EQ0, EQ1, EQ2, EQ3, EQ4, where EQ4 and EQ0 are reactant and product respectively) and 2 Transition states, TS0 (which connects EQ0 with EQ1) and TS1 which connects EQ2 with EQ3). But, we missed the other two TSs between EQ4–EQ3 and between EQ2–EQ1.
Which option would be, in principle, more effective to find a path? To do RePATH (in the same folder where DS-AFIR was performed) or to calculate new DS-AFIRs to find the missing TSs from the each pair of EQ structures?
Also, another thing that caught our attention in this particular DS-AFIR output was the highlighted with a # 'hashtag' word #RIDGE-POINT (we could find this in the manual) that appears instead of "#STEP N" in the main output XXX.log. This word seems to mark the final structure in the PATH LUP (after the structure marked with that word, the next line is this: "DS-AFIR COMPLETED, LUP-PATH OPTIMIZATION CARRIED OUT" ). We are wondering how we find the TS rather than the "ridge point"?
Thank you very much.
Re:
Satoshi Maeda
In general, I recommend to do RePATH in which the previous DS-AFIR results are read.
The ridge point will be used as the TS guess only when the DS = SingleSTEP option is used. Without this option, LUP is performed after the AFIR path is obtained, and energy maxima along the final LUP path are used as the TS guesses.