diff --git a/README.md b/README.md index 6156351..610ab20 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,12 @@ A unified interface for molecular harmonic vibrational frequency calculations. The UniMoVib program was originally written by Wenli Zou in FORTRAN 77 during 2014 and 2015 at Southern Methodist University (SMU), Dallas, Texas, within the framework of the LocalMode (now LModeA) program of the Computational and Theoretical Chemistry Group (CATCO) of SMU. This work was supported by the NSF grants CHE 1152357 and CHE 1464906. Guidance from the late Dr. Dieter Cremer is acknowledged. After being rewritten in Fortran 90 in the spring of 2017, UniMoVib has been released as a stand-alone program. ## Latest Versions +Version 1.5.0 (Jan/19/2023). + +1. Checkdata is printed for [BDF](http://bdf.theochem.cn:7226). +2. Gradient information is printed for [Gaussian](http://www.gaussian.com/), [CFour](http://www.cfour.de/) (analytical Hessian only), and UniMoVib. +3. The format of the UniMoVib data file has been updated. + Version 1.4.4 (Dec/28/2021). 1. Improved GCC version 10 (gfortran) compatibility. diff --git a/bin/win32.zip b/bin/win32.zip index f37963d..dd5b596 100644 Binary files a/bin/win32.zip and b/bin/win32.zip differ diff --git a/manual/manual-cn/manual-cn.pdf b/manual/manual-cn/manual-cn.pdf index 0147b75..9b26839 100644 Binary files a/manual/manual-cn/manual-cn.pdf and b/manual/manual-cn/manual-cn.pdf differ diff --git a/manual/manual-cn/manual-cn.tex b/manual/manual-cn/manual-cn.tex index caa178b..471a4af 100644 --- a/manual/manual-cn/manual-cn.tex +++ b/manual/manual-cn/manual-cn.tex @@ -8,7 +8,7 @@ \title{ \vspace{-3cm} \zihao{1}\hei\textsc{UniMoVib}使用手册 \\ -\vspace{1cm} \zihao{3}（1.4.4版）\vspace{20 mm} \\ +\vspace{1cm} \zihao{3}（1.5.0版）\vspace{20 mm} \\ \includegraphics[width=50.0mm]{fig/logo} \vspace{10 mm}} \author{ @@ -585,19 +585,20 @@ \chapter{附录} \label{part:appdx} \section{UniMoVib数据文件格式} \label{sec:almfmt} -\textsc{UniMoVib}数据文件是不区分大小写的纯文本自由格式，其中前三行的顺序是固定的，之后\verb|END|除外的数据区顺序任意。 +\textsc{UniMoVib}数据文件是不区分大小写的纯文本自由格式，其中第一行为标题行，之后关键词\verb|END|除外的数据区顺序任意。 +所有的关键词必须从第一列开始。注释行以空格或“!”打头，但不能出现在\{关键词$\ldots$数据\}之内。 -\noindent -1.1.0版(2020.10.24) +\vskip 0.5em \noindent +1.2.0版（2023.01.19） \begin{colorboxed}[oval=false,boxcolor=blue!75!black,bgcolor=blue!5!white] \ttfamily \begin{lstlisting} -（一个标题行） +（一个标题行；如果续行，则其余标题行的第一列必须是空格或“!”） NATM （一个正整数） AMASS - （NATM个原子质量；若不提供原子质量，可以忽略这部分内容或改用NOMASS， - 程序默认用最丰同位素质量） + （NATM个原子质量；若不提供原子质量，可以忽略这部分内容或改用NOMASS，程序将使用 + 最丰同位素质量） ZA （NATM个核电荷数） XYZ @@ -607,13 +608,13 @@ \section{UniMoVib数据文件格式} \label{sec:almfmt} FFX （3NATM*3NATM，Hessian矩阵元；若为下三角矩阵数据，改用 FFXLT） APT - （3*3NATM，APT数据值；若不提供 APT数据，可以忽略这部分内容或改用 NOAPT） + （3*3NATM，APT数据值用于计算红外强度；若不提供 APT数据，可以忽略这部分内容 + 或改用 NOAPT） DPR - （6*3NATM，极化率导数数据值；若以9*3NATM格式给出，改用DPRSQ； + （6*3NATM，极化率导数数据值用于计算拉曼活性；若以9*3NATM格式给出，改用DPRSQ； 若不提供 DPR数据，可以忽略这部分内容或改用 NODPR） END （可选的结束行） \end{lstlisting}\end{colorboxed} -目前\verb|GRD|数据仅用于某些外部程序。 \section{开发人员：到其它量子化学程序的接口} \label{sec:interface} @@ -627,29 +628,33 @@ \section{开发人员：到其它量子化学程序的接口} \label{sec:interfa ! Read NAtm from XXXX output !----------------------------------------------------------------------- subroutine RdNAtmXXXX(ifchk,NAtm,tag,ctmp) -implicit real(kind=8) (a-h,o-z) -character*100 :: ctmp -character*100 :: tag -(...) -return + implicit real(kind=8) (a-h,o-z) + character*100 :: ctmp + character*100 :: tag + (...) + return end \end{lstlisting}\end{colorboxed} \begin{enumerate} -\item[2] 读取原子单位的直角坐标（\verb|XYZ|），原子核电荷数（\verb|ZA|），原子单位的原子量或同位素质量（\verb|AMass|；可选），原子单位的直角坐标力常数矩阵（\verb|FFx|，必须是未做质量加权的对称方阵），原子单位的偶极矩梯度（即原子极化张量，\verb|APT|；可选），以及原子单位的极化率导数（\verb|DPol|；可选）。该子例程在\verb|subroutine RdData1|中调用。例如： +\item[2] 读取原子单位的直角坐标（\verb|XYZ|），原子核电荷数（\verb|ZA|），原子单位的原子量或同位素质量（\verb|AMass|；可选）， + 原子单位的直角梯度（\verb|Grad|；可选；读入后把\verb|IGrd|设为1），原子单位的直角坐标力常数矩阵（\verb|FFx|，必须是未做质量加权的对称方阵）， + 原子单位的偶极矩梯度（即原子极化张量，\verb|APT|；可选；读入后把\verb|Infred|设为1）， + 以及原子单位的极化率导数（\verb|DPol|；可选；读入后把\verb|IRaman|设为1）。该子例程在\verb|subroutine RdData1|中调用。例如： \end{enumerate} \begin{colorboxed}[oval=false,boxcolor=blue!75!black,bgcolor=blue!5!white] \begin{lstlisting}[language={[90]Fortran}] !----------------------------------------------------------------------- ! Read data from XXXX output !----------------------------------------------------------------------- -subroutine RdXXXX(ifchk,tag,ctmp,NAtm,AMass,ZA,XYZ,FFx,APT,DPol) -implicit real(kind=8) (a-h,o-z) -real(kind=8) :: AMass(NAtm),ZA(NAtm),XYZ(3,NAtm),FFx(3*NAtm,3*NAtm), & - APT(3,3*NAtm),DPol(6,3*NAtm) -character*100 :: ctmp -character*100 :: tag -(...) -return +subroutine RdXXXX(ifchk,tag,ctmp,IGrd,Infred,IRaman,NAtm,AMass,ZA,XYZ, & + Grad,FFx,APT,DPol) + implicit real(kind=8) (a-h,o-z) + real(kind=8) :: AMass(NAtm),ZA(NAtm),XYZ(3,NAtm),Grad(3,NAtm), & + FFx(3*NAtm,3*NAtm),APT(3,3*NAtm),DPol(6,3*NAtm) + character*100 :: ctmp + character*100 :: tag + (...) + return end \end{lstlisting}\end{colorboxed} \begin{enumerate} diff --git a/manual/manual-en/manual-en.pdf b/manual/manual-en/manual-en.pdf index 87fa7e4..36f4eba 100644 Binary files a/manual/manual-en/manual-en.pdf and b/manual/manual-en/manual-en.pdf differ diff --git a/manual/manual-en/manual-en.tex b/manual/manual-en/manual-en.tex index 903c1a8..03b1acc 100644 --- a/manual/manual-en/manual-en.tex +++ b/manual/manual-en/manual-en.tex @@ -52,7 +52,7 @@ \begin{document} \title{User's Guide of the Program \textsc{UniMoVib} \\ -\vspace{10 mm} (Ver. 1.4.4) \vspace{30 mm} \\ +\vspace{10 mm} (Ver. 1.5.0) \vspace{30 mm} \\ \includegraphics[width=50.0mm]{fig/logo} \vspace{30 mm} } \date{\today} @@ -838,11 +838,14 @@ \section{Appendix} \label{part:appdx} \subsection{Format of the UniMoVib data file} \label{sec:almfmt} \index{{UniMoVib format}@{UniMoVib format}} -The \textsc{UniMoVib} data file is a case-insensitive ASCII one in free format, where the order of the first three lines is fixed and the other data blocks except \verb|END| can be in any order. -\begin{Verbatim}[frame=single,label=Format(Ver.1.1.0 2020.10.24),labelposition=topline,rulecolor=\color{blue},fontsize=\small,baselinestretch=1.0] -(One title line) +The \textsc{UniMoVib} data file is a case-insensitive ASCII one in free format, where the first line is reserved for the title and the other data blocks except the keyword \verb|END| may be defined in any order. +All the keywords must start from the first column. The comment line starts with a space or !, but cannot be in the \{Keyword$\ldots$Data\} block. + +\begin{Verbatim}[frame=single,label=Format(Ver.1.2.0 2023.01.19),labelposition=topline,rulecolor=\color{blue},fontsize=\small,baselinestretch=1.0] +(One title line; for continuous multiple title lines, the first column must + be a space or !) NATM - (An positive integer) + (A positive integer) AMASS (NATM number of atomic masses; Use NOMASS instead or omit this block if no atomic masses provided; @@ -859,15 +862,14 @@ \subsection{Format of the UniMoVib data file} \label{sec:almfmt} (3NATM*3NATM elements of Hessian matrix; Use FFXLT instead for L.T. matrix) APT - (3*3NATM elements of APT data; + (3*3NATM elements of APT data for IR intensity; Use NOAPT instead or omit this block if no APT data provided) DPR - (6*3NATM elements of polarizability derivatives; + (6*3NATM elements of polarizability derivatives for Raman activities; Use DPRSQ instead if in the form of 9*3NATM; Use NODPR instead or omit this block if no DPR data provided) END (An optional end line) \end{Verbatim} -At present only some external programs need the \verb|GRD| data. \pagebreak{} @@ -883,30 +885,33 @@ \subsection{For developers: interface to other QC programs} \label{sec:interface ! Read NAtm from XXXX output !----------------------------------------------------------------------- subroutine RdNAtmXXXX(ifchk,NAtm,tag,ctmp) -implicit real(kind=8) (a-h,o-z) -character*100 :: ctmp -character*100 :: tag -(...) -return + implicit real(kind=8) (a-h,o-z) + character*100 :: ctmp + character*100 :: tag + (...) + return end \end{lstlisting} -\item Read in Cartesian coordinates in \emph{a.u.} (\verb|XYZ|), atomic nuclear charge number (\verb|ZA|), atomic or isotopic masses in \emph{a.u.} (\verb|AMass|; optional), Cartesian force constant matrix in \emph{a.u.} -(\verb|FFx|; a mass-unweighted square matrix), dipole moment gradient (\emph{i.e.} atomic polar tensor; APT) in \emph{a.u.} -(\verb|APT|; optional), and polarizability derivatives in \emph{a.u.} -(\verb|DPol|; optional). This subroutine is called in \verb|subroutine RdData1| . Example: +\item Read in Cartesian coordinates in \emph{a.u.} (\verb|XYZ|), atomic nuclear charge number (\verb|ZA|), atomic or isotopic masses in \emph{a.u.} + (\verb|AMass|; optional), Cartesian gradients in \emph{a.u.} (\verb|Grad|; optional; set \verb|IGrd=1| after reading the data), + Cartesian force constant matrix in \emph{a.u.} (\verb|FFx|; a mass-unweighted square matrix), + dipole moment gradient (\emph{i.e.} atomic polar tensor; APT) in \emph{a.u.} (\verb|APT|; optional; set \verb|Infred=1| after reading the data), + and polarizability derivatives in \emph{a.u.} (\verb|DPol|; optional; set \verb|IRaman=1| after reading the data). + This subroutine is called in \verb|subroutine RdData1| . Example: \begin{lstlisting}[language={[90]Fortran}] !----------------------------------------------------------------------- ! Read data from XXXX output !----------------------------------------------------------------------- -subroutine RdXXXX(ifchk,tag,ctmp,NAtm,AMass,ZA,XYZ,FFx,APT,DPol) -implicit real(kind=8) (a-h,o-z) -real(kind=8) :: AMass(NAtm),ZA(NAtm),XYZ(3,NAtm),FFx(3*NAtm,3*NAtm), & - APT(3,3*NAtm),DPol(6,3*NAtm) -character*100 :: ctmp -character*100 :: tag -(...) -return +subroutine RdXXXX(ifchk,tag,ctmp,IGrd,Infred,IRaman,NAtm,AMass,ZA,XYZ, & + Grad,FFx,APT,DPol) + implicit real(kind=8) (a-h,o-z) + real(kind=8) :: AMass(NAtm),ZA(NAtm),XYZ(3,NAtm),Grad(3,NAtm), & + FFx(3*NAtm,3*NAtm),APT(3,3*NAtm),DPol(6,3*NAtm) + character*100 :: ctmp + character*100 :: tag + (...) + return end \end{lstlisting} diff --git a/src/interface.f90 b/src/interface.f90 index 0716a08..eec0dd6 100644 --- a/src/interface.f90 +++ b/src/interface.f90 @@ -81,7 +81,7 @@ subroutine RdNAtm1(iinp,idt0,idt1,Intact,IOP,NAtm,tag,ctmp) return case(-2) ! UniMoVib (alm) - call RdNAtmALM(idt0,NAtm) + call RdNAtmALM(idt0,NAtm,ctmp) case(-3) ! xyz call RdNAtmXYZ(idt0,NAtm,ctmp) @@ -187,10 +187,10 @@ subroutine RdNAtm1(iinp,idt0,idt1,Intact,IOP,NAtm,tag,ctmp) ! IfFXX = .false. : do not read Hessian, and an approximate Hessian will be constructed later (for expert only!) ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,NAtm,tag,ctmp,AMass,ZA,XYZ,Grd,FFx,APT,DPol, & - Scr1,Scr2,Scr3,ScrD) +subroutine RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,ifbdfchk,NAtm,tag,ctmp,AMass,ZA,XYZ,Grd,FFx, & + APT,DPol,Scr1,Scr2,Scr3,ScrD) implicit real(kind=8) (a-h,o-z) -logical :: Intact, IfFXX +logical :: Intact, ifbdfchk, IfFXX dimension :: IOP(*) real(kind=8) :: AMass(*), ZA(*), XYZ(*), Grd(*), FFx(*), APT(*), DPol(*) real(kind=8) :: Scr1(*), Scr2(*), Scr3(*), ScrD(*) ! Scr1,2,3(NSS) and ScrD(2*NSS) @@ -200,6 +200,7 @@ subroutine RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,N Infred=0 IRaman=0 IGrd =0 +ifbdfchk = .false. IfFXX = IOP(10) == 0 select case(IOP(1)) @@ -208,7 +209,7 @@ subroutine RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,N return case(-2) ! UniMoVib (ALM); the size of ScrD should be 9*NAtm3 at least - call RdALMode(idt0,Intact,Infred,IRaman,IGrd,NAtm,ctmp,AMass,ZA,XYZ,Grd,FFx,APT,DPol,ScrD) + call RdALMode(idt0,Intact,Infred,IRaman,IGrd,ifbdfchk,NAtm,ctmp,AMass,ZA,XYZ,Grd,FFx,APT,DPol,ScrD) ! the most abundant isotopic masses are assumed if(AMass(1) < 0.d0) call MasLib(0,NAtm,AMass,ZA) @@ -234,7 +235,7 @@ subroutine RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,N call RdORCA(idt0,ctmp,Intact,Infred,IRaman,NAtm,AMass,ZA,XYZ,FFx,APT,DPol,Scr1) case(5) ! CFour - call RdCFour(idt0,idt1,tag,ctmp,Intact,IfFXX,Infred,NAtm,AMass,ZA,XYZ,FFx,APT,Scr1,Scr2,Scr3,ScrD) + call RdCFour(idt0,idt1,tag,ctmp,Intact,IfFXX,Infred,IGrd,NAtm,AMass,ZA,XYZ,Grd,FFx,APT,Scr1,Scr2,Scr3,ScrD) case(6) ! Molpro call RdMolp(idt0,tag,ctmp,Intact,Infred,NAtm,AMass,ZA,XYZ,FFx,APT) @@ -404,19 +405,33 @@ subroutine RdBMat(iout,ibmt,ctmp,Intact,NAtm3,IOP,FFx) ! Read NAtm from UniMoVib (ALM) data file ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine RdNAtmALM(ifchk,NAtm) -implicit real(kind=8) (a-h,o-z) -logical :: Intact +subroutine RdNAtmALM(ifchk,NAtm,ctmp) + implicit real(kind=8) (a-h,o-z) + character*100 :: ctmp -i=0 + i=0 -rewind(ifchk) -read(ifchk,"(/)",err=100,end=100) -read(ifchk,*,err=100,end=100) i -100 NAtm=i + rewind(ifchk) + ! the first title line + read(ifchk,*,err=100,end=100) + do while(.true.) + read(ifchk,"(a100)",err=100,end=100) ctmp + if(len_trim(ctmp) == 0) cycle + if(ctmp(1:1) == ' ' .or. ctmp(1:1) == '!') cycle -return -end + call charl2u(ctmp) + + ! read natm + if(index(ctmp,"NATM") == 1) then + read(ifchk,*,err=100,end=100) i + exit + end if + end do + + 100 NAtm=i + + return +end subroutine RdNAtmALM !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! @@ -1251,15 +1266,15 @@ subroutine RdNAtmACES(ifchk,Intact,NAtm,tag,ctmp) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! -! Read data from UniMoVib (ALM) data file (Ver. 1.1) +! Read data from UniMoVib (ALM) data file (Ver. 1.2) ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine RdALMode(ifchk,Intact,Infred,IRaman,IGrd,NAtm,ctmp,AMass,ZA,XYZ,Grd,FFx,APT,DPol,Scr) +subroutine RdALMode(ifchk,Intact,Infred,IRaman,IGrd,ifbdfchk,NAtm,ctmp,AMass,ZA,XYZ,Grd,FFx,APT,DPol,Scr) implicit real(kind=8) (a-h,o-z) parameter(ang2au=1.d0/0.52917720859d0) real(kind=8) :: AMass(NAtm),ZA(NAtm),XYZ(NAtm*3),Grd(NAtm*3),FFx(NAtm*NAtm*9),APT(NAtm*9),DPol(6,NAtm*3),Scr(*) character*100 :: ctmp - logical :: Intact + logical :: Intact, ifbdfchk NAtm3 = NAtm * 3 NAtm9 = NAtm * 9 @@ -1276,59 +1291,89 @@ subroutine RdALMode(ifchk,Intact,Infred,IRaman,IGrd,NAtm,ctmp,AMass,ZA,XYZ,Grd,F DPol = 0.d0 rewind(ifchk) - read(ifchk,"(//)",err=2010,end=2010) + + ! saved by BDFOPT? + read(ifchk,"(a100)",err=2000,end=2000) ctmp + call charl2u(ctmp) + ! UniMoVib data file generated by BDFOPT(Check) + if(index(ctmp,"UNIMOVIB") > 0 .and. index(ctmp,"BDFOPT(CHECK)") > 0) ifbdfchk = .true. + + read(ifchk,*,err=2010,end=2010) do while(.true.) read(ifchk,"(a100)",err=2020,end=1000) ctmp + if(len_trim(ctmp) == 0) cycle + if(ctmp(1:1) == ' ' .or. ctmp(1:1) == '!') cycle + call charl2u(ctmp) + ! read natm + if(index(ctmp,"NATM") == 1) then + read(ifchk,*) + ! read mass - if(index(ctmp,"MASS") > 0 .and. index(ctmp,"NOMASS") == 0) then + else if(index(ctmp,"MASS") == 1) then read(ifchk,*,err=2030,end=2030)(AMass(i),i=1,NAtm) + else if(index(ctmp,"NOMASS") == 1) then + AMass = 0.d0 + AMass(1)=-1.d0 + ! read IZ - else if(index(ctmp,"ZA") > 0 .and. index(ctmp,"XYZANG") == 0) then + else if(index(ctmp,"ZA") == 1) then read(ifchk,*,err=2040,end=2040)(ZA(i),i=1,NAtm) ! read Cartesian coordinates in a.u. (XYZ) or in Ang. (XYZANG) - else if(index(ctmp,"XYZ") > 0) then + else if(index(ctmp,"XYZ") == 1) then Scr(1)=1.d0 - if(index(ctmp,"XYZANG") > 0) Scr(1)=ang2au + if(index(ctmp,"XYZANG") == 1) Scr(1)=ang2au read(ifchk,*,err=2050,end=2050)(XYZ(i),i=1,NAtm3) ! Ang --> a.u. call AScale(NAtm3,Scr(1),XYZ,XYZ) + ! read GRD in a.u. + else if(index(ctmp,"GRD") == 1) then + IGrd = 1 + read(ifchk,*,err=2090,end=2090)(Grd(i),i=1,NAtm3) + + else if(index(ctmp,"NOGRD") == 1) then + IGrd = 0 + Grd = 0.d0 + ! read square or L.T. FFX in a.u. - else if(index(ctmp,"FFX") > 0) then - if(index(ctmp,"FFXLT") == 0) then - read(ifchk,*,err=2060,end=2060)(FFX(i),i=1,NSS) - else + else if(index(ctmp,"FFX") == 1) then + if(index(ctmp,"FFXLT") == 1) then read(ifchk,*,err=2060,end=2060)(Scr(i),i=1,NTT) call LT2Sqr(NAtm3,Scr,FFx) + else + read(ifchk,*,err=2060,end=2060)(FFX(i),i=1,NSS) end if ! read APT in a.u. - else if(index(ctmp,"APT") > 0 .and. index(ctmp,"NOAPT") == 0) then + else if(index(ctmp,"APT") == 1) then Infred = 1 read(ifchk,*,err=2070,end=2070)(APT(i),i=1,NAtm9) + else if(index(ctmp,"NOAPT") == 1) then + Infred = 0 + APT = 0.d0 + ! read DPR in a.u. - else if(index(ctmp,"DPR") > 0 .and. index(ctmp,"NODPR") == 0) then + else if(index(ctmp,"DPR") == 1) then IRaman = 1 - if(index(ctmp,"DPRSQ") == 0) then - read(ifchk,*,err=2080,end=2080)((DPol(j,i),j=1,6),i=1,NAtm3) - else + if(index(ctmp,"DPRSQ") == 1) then read(ifchk,*,err=2080,end=2080)(Scr(i),i=1,NAtm3*9) call S9to6(NAtm3,Scr,DPol) + else + read(ifchk,*,err=2080,end=2080)((DPol(j,i),j=1,6),i=1,NAtm3) end if - ! read GRD in a.u. - else if(index(ctmp,"GRD") > 0 .and. index(ctmp,"NOGRD") == 0) then - IGrd = 1 - read(ifchk,*,err=2090,end=2090)(Grd(i),i=1,NAtm3) + else if(index(ctmp,"NODPR") == 1) then + IRaman = 0 + DPol = 0.d0 ! read END - else if(index(ctmp,"END") > 0) then + else if(index(ctmp,"END") == 1) then exit end if @@ -1336,6 +1381,7 @@ subroutine RdALMode(ifchk,Intact,Infred,IRaman,IGrd,NAtm,ctmp,AMass,ZA,XYZ,Grd,F 1000 return + 2000 call XError(Intact,"Please check TITLE of the UniMoVib data file!") 2010 call XError(Intact,"Please check NATM data!") 2020 call XError(Intact,"Please check UniMoVib data file!") 2030 call XError(Intact,"Please check MASS data!") @@ -1345,7 +1391,7 @@ subroutine RdALMode(ifchk,Intact,Infred,IRaman,IGrd,NAtm,ctmp,AMass,ZA,XYZ,Grd,F 2070 call XError(Intact,"Please check APT data!") 2080 call XError(Intact,"Please check DPR/DPRSQ data!") 2090 call XError(Intact,"Please check GRD data!") -end +end subroutine RdALMode !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! @@ -1472,6 +1518,8 @@ subroutine RdGauss(ifchk,tag,ctmp,Intact,IRdMas,Infred,IRaman,IGrd,NAtm,AMass,ZA 161 read(ifchk,"(a100)",end=200)ctmp if(index(ctmp,tag)==0)goto 161 read(ifchk,"(5e16.8)")(Grd(i),i=1,NAtm3) +! in fact they are forces (negative of the gradients) +Grd(1:NAtm3) = -Grd(1:NAtm3) IGrd= 1 ! read FFx (a.u.) @@ -1749,9 +1797,9 @@ logical function ChkAFRQCfour(ifchk,tag,ctmp,Intact) ! Read data from CFour *.out ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine RdCFour(ifchk,igeom,tag,ctmp,Intact,IfFXX,Infred,NAtm,AMass,ZA,XYZ,FFx,APT,Scr1,Scr2,Scr3,Scr4) +subroutine RdCFour(ifchk,igeom,tag,ctmp,Intact,IfFXX,Infred,IGrd,NAtm,AMass,ZA,XYZ,Grd,FFx,APT,Scr1,Scr2,Scr3,Scr4) implicit real(kind=8) (a-h,o-z) - real(kind=8) :: AMass(*),ZA(*),XYZ(3,*),FFx(NAtm*3,*),APT(3,*),Scr1(*),Scr2(*),Scr3(*),Scr4(*) + real(kind=8) :: AMass(*),ZA(*),XYZ(3,*),Grd(3,*),FFx(NAtm*3,*),APT(3,*),Scr1(*),Scr2(*),Scr3(*),Scr4(*) character*100 :: tag,ctmp logical :: Intact,IfFXX,itrue,ChkAFRQCfour @@ -1761,7 +1809,7 @@ subroutine RdCFour(ifchk,igeom,tag,ctmp,Intact,IfFXX,Infred,NAtm,AMass,ZA,XYZ,FF itrue = itrue .and. ChkAFRQCfour(ifchk,tag,ctmp,Intact) if(itrue) then - call RdCFourA(ifchk,igeom,tag,ctmp,Intact,Infred,NAtm,AMass,ZA,XYZ,FFx,APT,Scr1) + call RdCFourA(ifchk,igeom,tag,ctmp,Intact,Infred,IGrd,NAtm,AMass,ZA,XYZ,Grd,FFx,APT,Scr1) else Infred= 0 call RdCFourN(ifchk,tag,ctmp,Intact,IfFXX,NAtm,AMass,ZA,XYZ,FFx,Scr1,Scr2,Scr3,Scr4) @@ -1915,9 +1963,9 @@ subroutine RdCFourN(ifchk,tag,ctmp,Intact,IfFXX,NAtm,AMass,ZA,XYZ,FFx,SC1,SC2,SC ! Read data from CFour *.out & GRD. For analytic frequency only. ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine RdCFourA(ifchk,igeom,tag,ctmp,Intact,Infred,NAtm,AMass,ZA,XYZ,FFx,APT,Scr) +subroutine RdCFourA(ifchk,igeom,tag,ctmp,Intact,Infred,IGrd,NAtm,AMass,ZA,XYZ,Grd,FFx,APT,Scr) implicit real(kind=8) (a-h,o-z) -real(kind=8) :: AMass(*),ZA(*),XYZ(3,*),FFx(NAtm*3,*),APT(3,*),Scr(*) +real(kind=8) :: AMass(*),ZA(*),XYZ(3,*),Grd(3,*),FFx(NAtm*3,*),APT(3,*),Scr(*) character*100 :: ctmp character*77 :: tag,tag1 logical :: Intact,ifc1 @@ -1934,6 +1982,12 @@ subroutine RdCFourA(ifchk,igeom,tag,ctmp,Intact,Infred,NAtm,AMass,ZA,XYZ,FFx,APT read(igeom,*)ZA(i),(XYZ(j,i),j=1,3) end do call ACopy(NAtm3,XYZ,Scr) +! read Cartesian gradients +do i=1,NAtm + read(igeom,*) x,(Grd(j,i),j=1,3) + if(nint(x) /= nint(ZA(i))) call XError(Intact,"In the GRD file, ZA1 .ne. ZA2!") +end do +IGrd=1 ! C1 symmetry? tag=' The computational point group is ' @@ -4164,7 +4218,7 @@ subroutine SavALM(ifalm,ifloc,ialm,iloc,Infred,IRaman,IGrad,NAtm,NAtm3,NVib,AMas if(ifalm) then rewind(ialm) - write(ialm,"(' UniMoVib DATA FILE (THIS TITLE CAN BE MODIFIED)')") + write(ialm,"(' UniMoVib DATA FILE',/,' (THIS TITLE CAN BE MODIFIED)')") write(ialm,"('NATM')") write(ialm,"(i5)")NAtm @@ -4362,7 +4416,7 @@ subroutine SavGJF(igau,NAtm,NAtm3,NVib,ifreq,Rslt,XYZ,ZA,AMass,FFx,AL,Scr1,Scr2, write(igau,"('%mem=2GB')") write(igau,"('%nprocshared=2')") - write(igau,"('#p b3lyp/genecp iop(1/10=1) opt(ts,ReadFC,cartesian,noeigentest)')") + write(igau,"('#p b3lyp/genecp iop(1/10=1) opt(ts,ReadFC,cartesian,noeigentest,nomicro)')") write(igau,"(/,'This is a templet for TS optimization. Do not use it for ONIOM!',/)") totel=sum(ZA) diff --git a/src/main.f90 b/src/main.f90 index 985b065..e358e06 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -24,17 +24,17 @@ program UniMoVib character*5 :: ver character*12 :: dat character*200 :: ctmp, cname, tag - logical :: Intact, ifopen + logical :: Intact, ifopen, ifbdfchk allocatable :: AMass(:), ZA(:), XYZ(:), Grd(:), FFx(:), APT(:), DPol(:), AL(:), Rslt(:) allocatable :: Scr1(:), Scr2(:), Scr3(:), Scr4(:), Work(:) logical, allocatable :: LScr(:) integer, allocatable :: IScr(:) character*1, allocatable :: CScr(:) integer, allocatable :: subsystem_idx(:),flags(:) - allocatable :: AMass_sub(:), XYZ_sub(:), ZA_sub(:) + allocatable :: AMass_sub(:), XYZ_sub(:), ZA_sub(:) - ver="1.4.4" - dat="Dec 28, 2021" + ver="1.5.0" + dat="Jan 19, 2023" !----------------------------------------------------------------------- ! 1. Assign I/O @@ -98,8 +98,8 @@ program UniMoVib if(ierr /= 0) call XError(Intact,"Insufficient Memory (3)!") APT=0.d0 - call RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,NAtm,tag,ctmp,AMass,ZA,XYZ,Grd,FFx,APT,DPol, & - Scr1,Scr2,Scr3,Work) + call RdData1(iinp,iout,idt0,idt1,idt2,ibmt,Intact,IOP,Infred,IRaman,IGrd,ifbdfchk,NAtm,tag,ctmp,AMass,ZA,XYZ,Grd,FFx, & + APT,DPol,Scr1,Scr2,Scr3,Work) ! read atomic masses from input call RdIsot(iinp,iout,Intact,IOP(4),NAtm,ctmp,AMass) @@ -115,8 +115,8 @@ program UniMoVib ! 7. Solve Secular equation in Cartesian coordinates ! Symmetry is also analyzed therein. !----------------------------------------------------------------------- - call SolvSec(iinp,iout,idt0,irep,ireo,iudt,imdn,iloc,igau,imdf,Intact,IOP,Infred,IRaman,IGrd,NAtm,NVib,ctmp,AMass,ZA,XYZ, & - Grd,FFx,APT,DPol,AL,Rslt,LScr,IScr,CScr,Scr1,Scr2,Scr3,Scr4,Work) + call SolvSec(iinp,iout,idt0,irep,ireo,iudt,imdn,iloc,igau,imdf,Intact,IOP,Infred,IRaman,IGrd,ifbdfchk,NAtm,NVib,ctmp, & + AMass,ZA,XYZ,Grd,FFx,APT,DPol,AL,Rslt,LScr,IScr,CScr,Scr1,Scr2,Scr3,Scr4,Work) !----------------------------------------------------------------------- ! 8. GSVA - Generalized Subsystem Vibrational Analysis @@ -137,6 +137,11 @@ program UniMoVib call GSVA_engine(iout,isva,imds,IOP,NAtm,subsystem_idx,NAtm_sub,AMass_sub,XYZ_sub,ZA_sub,FFx) end if +!----------------------------------------------------------------------- +! 9. print gradient information +!----------------------------------------------------------------------- + if(IGrd /= 0) call GradInfo(iout,ifbdfchk,NAtm,NAtm3,NVib,ZA,Grd,FFx,AL,Rslt,Scr1,Scr2,ctmp) + !----------------------------------------------------------------------- ! 99. the last step !----------------------------------------------------------------------- diff --git a/src/rw.f90 b/src/rw.f90 index 71acb06..95047ca 100644 --- a/src/rw.f90 +++ b/src/rw.f90 @@ -772,5 +772,52 @@ subroutine dumpir(ida1,ida2,c4) return end +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! +! print convergence information of gradients. +! +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +subroutine prtconv(iout,grslt) + implicit real(kind=8) (a-h,o-z) + parameter(tolx=2.5d-3, tolg=5.0d-4, tole=5.0d-6) + dimension :: grslt(5) + + write(iout,"(/,' Convergence of gradients',/,34x,'Value Tolerance Converged?')") + nyes = 0 + if(grslt(1) <= tolx*1.6d0) then + write(iout,"(' Maximum Delta-X',8x,2f14.6,12x,'Yes')") grslt(1), tolx*1.6d0 + nyes = nyes + 1 + else + write(iout,"(' Maximum Delta-X',8x,2f14.6,13x,'No')") grslt(1), tolx*1.6d0 + end if + if(grslt(2) <= tolx) then + write(iout,"(' RMS Delta-X',8x,2f14.6,12x,'Yes')") grslt(2), tolx + nyes = nyes + 1 + else + write(iout,"(' RMS Delta-X',8x,2f14.6,13x,'No')") grslt(2), tolx + end if + if(grslt(3) <= tolg*1.6d0) then + write(iout,"(' Maximum Force',8x,2f14.6,12x,'Yes')") grslt(3), tolg*1.6d0 + nyes = nyes + 1 + else + write(iout,"(' Maximum Force',8x,2f14.6,13x,'No')") grslt(3), tolg*1.6d0 + end if + if(grslt(4) <= tolg) then + write(iout,"(' RMS Force',8x,2f14.6,12x,'Yes')") grslt(4), tolg + nyes = nyes + 1 + else + write(iout,"(' RMS Force',8x,2f14.6,13x,'No')") grslt(4), tolg + end if + if(grslt(5) <= tole) then + write(iout,"(' Expected Delta-E',8x,2d14.2,12x,'Yes')") grslt(5), tole + nyes = nyes + 1 + else + write(iout,"(' Expected Delta-E',8x,2d14.2,13x,'No')") grslt(5), tole + end if + if(nyes /= 5) write(iout,"(/,' This geometry is not a stationary point.')") + + return +end subroutine prtconv + !--- END diff --git a/src/thermo.f90 b/src/thermo.f90 index 198c967..7c06b51 100644 --- a/src/thermo.f90 +++ b/src/thermo.f90 @@ -1,6 +1,7 @@ !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! Thermochemistry calculation. Input: +! ifbdfchk: print check data for bdf ! NAtom : #atoms ! NAtm3 : NAtom*3 ! NVib : #vib. modes @@ -12,7 +13,7 @@ ! ctmp : scratch for characters. ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine Thermochem(iinp,iout,Intact,NAtom,NAtm3,NVib,IFAtom,IExpt,PGNAME,AMass,XYZ,Freq,Sc1,Sc2,ctmp) +subroutine Thermochem(iinp,iout,Intact,ifbdfchk,NAtom,NAtm3,NVib,IFAtom,IExpt,PGNAME,AMass,XYZ,Freq,Sc1,Sc2,ctmp) implicit real(kind=8) (a-h,o-z) parameter(tolr=0.01d0) parameter( pi = acos(-1.0d0), & @@ -38,7 +39,7 @@ subroutine Thermochem(iinp,iout,Intact,NAtom,NAtm3,NVib,IFAtom,IExpt,PGNAME,AMas character*4 :: PGNAME(2),PG character*1 :: L2U real(kind=8) :: energy(3),entropy(4) -logical :: Intact,IFAtom,IFLin,IFRdT,IfRdP +logical :: Intact,ifbdfchk,IFAtom,IFLin,IFRdT,IfRdP ! Eel : Total electronic energy from Q.C. calculation (a.u.) ! temp : temperature (K) @@ -248,6 +249,13 @@ subroutine Thermochem(iinp,iout,Intact,NAtom,NAtm3,NVib,IFAtom,IExpt,PGNAME,AMas ' Sum of electronic and thermal Free Energies:',f20.6)")zpe+Eel,eu+Eel,eh+Eel,eg+Eel write(iout,"(1x,84('=') )") + ! print check data for bdf + if(ifbdfchk) then + call bdfchk(iout,.false.) + write(iout,"(' CHECKDATA:BDFOPT:THERMO: ',4f16.2)") zpe*H2kcal,eu*H2kcal,eh*H2kcal,eg*H2kcal + call bdfchk(iout,.true.) + end if + end do 1000 continue @@ -255,7 +263,7 @@ subroutine Thermochem(iinp,iout,Intact,NAtom,NAtm3,NVib,IFAtom,IExpt,PGNAME,AMas return 2000 call XError(Intact,"Please check your input of $Thermo!") -end +end subroutine Thermochem !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! diff --git a/src/util.f90 b/src/util.f90 index 0c1e2d7..2b74583 100644 --- a/src/util.f90 +++ b/src/util.f90 @@ -118,7 +118,7 @@ function EleMas(Mode,IZ) 266.119830000d0, 267.121790000d0, 268.125670000d0, 269.128630000d0, 270.133360000d0, 277.151900000d0,& 278.156310000d0, 281.164510000d0, 283.170540000d0, 285.177120000d0, 286.182210000d0, 289.190420000d0,& 290.195980000d0, 293.204490000d0, 294.210460000d0, 294.213920000d0, 295.000000000d0, 299.000000000d0/ - ! averaged isotopic masses + ! average isotopic masses data (amas2(i),i=1,maxza) / & 1.007940000d0, 4.002602000d0, 6.941000000d0, 9.012183100d0, 10.811000000d0, 12.010700000d0,& 14.006700000d0, 15.999400000d0, 18.998403163d0, 20.179700000d0, 22.989769280d0, 24.305000000d0,& @@ -1516,4 +1516,93 @@ subroutine CountIrrep(Intact,iport,NAtm3,NVib,NClass,Irreps,ModMap) return end +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! +! gradient information +! +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +subroutine GradInfo(iout,ifbdfchk,NAtm,NAtm3,NVib,za,grd,ffx,al,fcon,dx,scr,Elm) + implicit real(kind=8) (a-h,o-z) + logical :: ifbdfchk + dimension :: za(NAtm), grd(3,NAtm), ffx(*), al(NAtm3,*), fcon(NAtm3), dx(3,*), scr(NAtm3,NAtm3) + character*3 :: Elm + allocatable :: grslt(:) + + write(iout,"(//,1x,30('*'),/,' *** Gradient Information ***',/,1x,30('*'))") + + write(iout,"(/,' Cartesian gradients (a.u.)',/,1x,68('-'),/,3x, & + 'No. Atom ZA X Y Z',/,1x,68('-'))") + do i=1,NAtm + j = nint(za(i)) + call ElemZA(1,Elm,j) + write(iout,"(i6,4x,a3,1x,i5,8x,3f14.8)") i,Elm,j,grd(:,i) + end do + write(iout,"(1x,68('-'))") + + ! NewtonCRaphson step dX = -F^-1 * g, + ! see Eq. (3) in JCP, 111, 10806 (1999), where H --> F and f --> -g. + ! It is worse than the RFO (rational function optimization) step (see Eq. (9)) but is much simpler. + + ! F^-1 calculation: + ! From F L = M L E, we have FF = LL E LL' where FF = X F X, LL = X^-1 L, and X = M^-1/2. + ! Then F^-1 = L E^-1 L'. + ! However AL saves renormalized L0 = L mu^-1/2, and E^-1 = mu k^-1, so F^-1 = L0 k^-1 L0'. + do i = 1, NVib + x = sign(max(abs(fcon(i)), 1.0d-8), fcon(i)) + scr(:,i) = (1.0d0/x) * al(:,i) + end do + call DGEMM('N','T',NAtm3,NAtm3,NVib,1.d0,scr,NAtm3,al,NAtm3,0.d0,dx,NAtm3) ! F^-1 --> dx + call MMpyMF(NAtm3,NAtm3,1,dx,grd,scr) + call AScale(NAtm3,-1.0d0,scr,dx) + + allocate(grslt(5)) + grslt = 0.0d0 + + do i = 1, NAtm + do j = 1, 3 + grslt(1) = max(grslt(1),abs(dx(j,i))) + grslt(2) = grslt(2) + dx(j,i) * dx(j,i) + grslt(3) = max(grslt(3),abs(grd(j,i))) + grslt(4) = grslt(4) + grd(j,i) * grd(j,i) + end do + end do + grslt(2) = sqrt(grslt(2)/dble(NAtm3)) + grslt(4) = sqrt(grslt(4)/dble(NAtm3)) + ! dE: see Eq. (1) in JCP, 111, 10806 (1999), where f = -g. + call MMpyMF(1,NAtm3,NAtm3,dx,ffx,scr) + grslt(5) = abs(0.5d0*dotx(NAtm3,scr,dx) + dotx(NAtm3,grd,dx)) + + call prtconv(iout,grslt) + + ! print check data for bdf + if(ifbdfchk) then + call bdfchk(iout,.false.) + write(iout,"(' CHECKDATA:BDFOPT:CONVERGE:',4f16.6)") grslt(1:4) + write(iout,"(' CHECKDATA:BDFOPT:CONVERGE:',d16.2)") grslt(5) + call bdfchk(iout,.true.) + end if + + deallocate(grslt) + + return +end subroutine GradInfo + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! +! print BDF CHECKDATA block. +! +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +subroutine bdfchk(iout,ifend) + implicit real(kind=8) (a-h,o-z) + logical :: ifend + + if(.not. ifend) then + write(iout,"(/,1x,39('+'),' DATA CHECK ',40('+'))") + else + write(iout,"(1x,37('+'),' END DATA CHECK ',38('+')),/") + end if + + return +end subroutine bdfchk + !--- END diff --git a/src/vibfreq.f90 b/src/vibfreq.f90 index 4d5c6eb..8f4f5c9 100644 --- a/src/vibfreq.f90 +++ b/src/vibfreq.f90 @@ -152,15 +152,15 @@ subroutine GSVA_engine(iout,isva,imds,IOP,NAtm,subsystem_idx,NAtm_sub,AMass_sub, ! Solve Secular equation in Cartesian coordinates ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine SolvSec(iinp,iout,idt0,irep,ireo,iudt,imdn,iloc,igau,imdf,Intact,IOP,Infred,IRaman,IGrd,NAtm,NVib,ctmp,AMass, & - ZA,XYZ,Grd,FFx,APT,DPol,AL,Rslt,LScr,IScr,Irname,Scr1,Scr2,Scr3,Scr4,Work) +subroutine SolvSec(iinp,iout,idt0,irep,ireo,iudt,imdn,iloc,igau,imdf,Intact,IOP,Infred,IRaman,IGrd,ifbdfchk,NAtm,NVib,ctmp, & + AMass,ZA,XYZ,Grd,FFx,APT,DPol,AL,Rslt,LScr,IScr,Irname,Scr1,Scr2,Scr3,Scr4,Work) implicit real(kind=8) (a-h,o-z) dimension :: IOP(*),AMass(*),ZA(*),XYZ(*),Grd(*),FFx(*),APT(*),DPol(*),AL(*),Rslt(*),IScr(*),Scr1(*),Scr2(*) ! NOTE. Scr3, Scr4, and Work are not available if IOP(9) = 1. dimension :: Scr3(*),Scr4(*),Work(*) character*4 :: PGNAME(2),Irname(NAtm*3) character*100 :: ctmp - logical :: Intact,IFAtom,IfSimu,LScr(*) + logical :: Intact,ifbdfchk,IFAtom,IfSimu,LScr(*) IFAtom = IOP(1) == -1 NAtm3=3*NAtm @@ -211,7 +211,7 @@ subroutine SolvSec(iinp,iout,idt0,irep,ireo,iudt,imdn,iloc,igau,imdf,Intact,IOP, ! print normal modes results saved in Rslt IfSimu = IOP(1) == -3 .or. IOP(1) == -4 - call PrtNFq(iout,irep,imdf,Infred,IRaman,NAtm,NAtm3,NVib,IfSimu,IOP(2),IOP(10),IOP(16),ZA,AL,Rslt,Irname) + call PrtNFq(iout,irep,imdf,Infred,IRaman,ifbdfchk,NAtm,NAtm3,NVib,IfSimu,IOP(2),IOP(10),IOP(16),ZA,AL,Rslt,Irname) ! experimental frequencies if(IOP(3) == 1) call RdExFq(iinp,iout,irep,Intact,NAtm3,NVib,Rslt,Scr2,Irname,LScr,Work,ctmp) @@ -231,10 +231,10 @@ subroutine SolvSec(iinp,iout,idt0,irep,ireo,iudt,imdn,iloc,igau,imdf,Intact,IOP, 1000 continue ! Thermochemistry calculation. Frequencies are saved in Rslt(1:NVib,3/5) in a.u. - call Thermochem(iinp,iout,Intact,NAtm,NAtm3,NVib,IFAtom,IOP(3),PGNAME,AMass,XYZ,Rslt,Scr1,Scr2,ctmp) + call Thermochem(iinp,iout,Intact,ifbdfchk,NAtm,NAtm3,NVib,IFAtom,IOP(3),PGNAME,AMass,XYZ,Rslt,Scr1,Scr2,ctmp) return -end +end subroutine SolvSec !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! @@ -832,12 +832,12 @@ subroutine PrtNFq_gsva(iout,NAtm,NAtm3,NVib,IPrint,ZA,subsystem_idx,AL,Reslt,Ipy ! print results of normal modes ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -subroutine PrtNFq(iout,irep,imdf,Infred,IRaman,NAtm,NAtm3,NVib,IfSim,IPrint,IApprox,Ipyvibms,ZA,AL,Reslt,IRNAME) +subroutine PrtNFq(iout,irep,imdf,Infred,IRaman,ifbdfchk,NAtm,NAtm3,NVib,IfSim,IPrint,IApprox,Ipyvibms,ZA,AL,Reslt,IRNAME) implicit real(kind=8) (a-h,o-z) parameter(au2wn=5140.48714376d0,au2dy=15.56893d0,cf=31.22307d0,au2ang4=0.52917720859d0**4) real(kind=8) :: ZA(*),AL(3,NAtm,*),Reslt(NAtm3,*) character*4 :: IRNAME(NAtm3) - logical :: IfSim + logical :: ifbdfchk,IfSim ! read irreps rewind(irep) @@ -965,7 +965,17 @@ subroutine PrtNFq(iout,irep,imdf,Infred,IRaman,NAtm,NAtm3,NVib,IfSim,IPrint,IApp !write(imdf,"(i5,3x,i5)") NAtm,NVib end if + ! print check data for bdf + if(ifbdfchk) then + call bdfchk(iout,.false.) + write(iout,"(' CHECKDATA:BDFOPT:NVIB:',11x,i4)") Nvib + do i = 1, Nvib + write(iout,"(' CHECKDATA:BDFOPT:FREQ:',3f17.1,f17.2)") Reslt(i,3)*au2wn,Reslt(i,4)*cf*cf,Reslt(i,7)*au2ang4,Reslt(i,8) + end do + call bdfchk(iout,.true.) + end if + return -end +end subroutine PrtNFq !--- END