comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 22:25:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

HvD: Finally fixed the bug that caused the geometry optimizer to fail in

Browse source 
case of symmetry and constrained geometric parameters. The problem was that
if a symmetric starting structure was provided but parameters (bond lengths,
bond angles, or torsion angles) were declared constant so that a geometry
optimization would lower the point group the code failed. In those cases the
code would enforce the point group symmetry of the input structure and ignore
the constraints rather than adjust the point group symmetry to be compatible
with the constraints.

The changes implemented make the symmetry code aware of any constraints imposed
on the geometry. This is done by generating a numeric label for each atom that
is different for each type of and number of constraints that an atom appears in.
For two atoms to be symmetry equivalent we insist that they must:
- be at positions that are related by the symmetry operators
- have the same tag
- have the same numeric constraint label.
The 4 calculations in the test case should all behave in the same way and give
the same result.
This commit is contained in:
Huub Van Dam 2014-10-07 21:24:59 +00:00
parent 82329626c2
commit 04ab47402e
8 changed files with 12458 additions and 9829 deletions
Download patch file Download diff file Expand all files Collapse all files

37
QA/tests/h2o_opt/README
View file

@ -1,25 +1,30 @@
Constraint optimization
=======================
This input runs 3 calculations on the water molecule. In all 3 calculations
This input runs 4 calculations on the water molecule. In all 4 calculations
one OH bond length is constraint, the other one is optimized as well as the
bond angle. However the geometries are specified slightly differently:
h2o_c1: The symmetry option is used to turn all symmetry off.
The resulting calculation works as it should.
h2o_c1: The symmetry option is used to turn all symmetry off.
The resulting calculation works as it should.
h2o_c2v: The input geometry is symmetric but the constraint is asymmetric
and autosym is not turned off.
The resulting calculation optimizes everything even the constraint
bondlength. THIS IS WRONG!!!
h2o_c2v: The input geometry is symmetric but the constraint is asymmetric
and autosym is not turned off.
The code now correctly finds Cs symmetry and completes the
optimization is performed correctly.
h2o_cs: The input geometry is asymmetric otherwise everything is the same
as with the h2o_c2v geometry.
The resulting calculation correctly sees the two bondlengths as
different and applies the constraints in the correct way.
h2o_cs: The input geometry is asymmetric otherwise everything is the same
as with the h2o_c2v geometry.
The resulting calculation correctly sees the two bondlengths as
different and applies the constraints in the correct way.
If correct then all three calculations should arrive at the same answer.
At the moment they do not. The baseline in this directory has been manually
constructed by copying the geometry optimization output for the h2o_c1
geometry to replace the incorrect results for the h2o_c2v geometry. At the
moment this is the best guess for what the correct output should look like.
h2o_cart: The input geometry is specified using Cartesian coordinates and
the constraints are imposed using the geometry adjust; zcoord;
mechanism.
The code detects Cs symmetry correctly also in this case and
completes the geometry optimization correctly.
If correct then all four calculations should arrive at the same answer,
following the same optimization path. Making the auto-symmetry code aware
of any constraints on the geometry fixed the optimization problems we had
before.

45
QA/tests/h2o_opt/h2o_opt.nw
View file

@ -1,13 +1,6 @@
echo
start h2o_opt_dat
# NOTE
# ====
#
# This test case currently fails, but it should actually complete
# successfully or stop with a helpful error message. See the README
# file for more details.
geometry h2o_c1 units angstroms
symmetry c1
zmatrix
@ -38,25 +31,57 @@ end
geometry h2o_cs units angstroms
zmatrix
O
H 1 OH1
H 1 OH2 2 HOH
H 1 OH1
H1 1 OH2 2 HOH
variables
HOH 108.0
OH1 1.1
OH1 1.2
constants
OH2 1.2
end
end
geometry h2o_cart units angstrom
O 0.00000000 0.00000000 -0.14106846
H -0.97082039 0.00000000 0.56427384
H 0.97082039 0.00000000 0.56427384
end
geometry h2o_cart adjust
zcoord
bond 1 3 1.2 r constant
end
end
basis
* library 6-31G*
end
driver
clear
end
set geometry h2o_c1
task scf optimize
driver
clear
end
set geometry h2o_c2v
task scf optimize
driver
clear
end
set geometry h2o_cs
task scf optimize
driver
clear
end
set geometry h2o_cart
task scf optimize

21897
QA/tests/h2o_opt/h2o_opt.out
View file

File diff suppressed because it is too large Load diff

90
src/geom/geom.F
View file

@ -518,7 +518,7 @@ c
end
c
C> \brief Load a geometry from the RTDB
c
C>
C> Load a geometry from the RTDB with a specified name.
C> The name is used also to define the name of the geometry.
C>
@ -609,6 +609,10 @@ c
tmp(k:) = ':masses'
s = s .and. rtdb_get(rtdb, tmp, mt_dbl, max_cent,
$ geom_mass(1,geom))
tmp(k:) = ' '
tmp(k:) = ':atomct'
s = s .and. rtdb_get(rtdb, tmp, mt_dbl, max_cent,
$ geom_atomct(1,geom))
C new
tmp(k:) = ' '
tmp(k:) = ':inv nuc expon'
@ -1104,6 +1108,10 @@ c
tmp(k:) = ':masses'
s = s .and. rtdb_put(rtdb, tmp, mt_dbl, ncenter(geom),
$ geom_mass(1,geom))
tmp(k:) = ' '
tmp(k:) = ':atomct'
s = s .and. rtdb_put(rtdb, tmp, mt_dbl, ncenter(geom),
$ geom_atomct(1,geom))
C new
tmp(k:) = ' '
tmp(k:) = ':inv nuc expon'
@ -4478,11 +4486,11 @@ c
geom_atn_to_default_mass = .true.
c
end
c
C>
C> \brief Define the atomic masses of the centers in a geometry instance
c
C>
C> \return Return .true. if successfull, and .false. otherwise
c
C>
logical function geom_masses_set(geom, ncent, masses)
implicit none
#include "nwc_const.fh"
@ -4612,6 +4620,80 @@ c
mass = geom_mass(icent,geom)
c
end
C>
C> \brief Define the atom constraint type of the centers in a geometry
C> instance
C>
C> \return Return .true. if successfull, and .false. otherwise
C>
logical function geom_atomct_set(geom, ncent, atomct)
implicit none
#include "nwc_const.fh"
#include "geomP.fh"
#include "stdio.fh"
c
integer geom !< [Input] the geometry handle
integer ncent !< [Input] the number of centers
double precision atomct(ncent) !< [Input] the atom constraint type
!< on each center
c
integer i
c
logical geom_check_handle
external geom_check_handle
c
geom_atomct_set = geom_check_handle(geom, 'geom_atomct_set')
if (.not. geom_atomct_set) return
c
if (ncent.le.0) then
write(LuOut,*) ' geom_atomct_set: too few centers ',ncent,
$ names(geom)(1:lenn(geom))
geom_atomct_set = .false.
return
else if (ncent.gt.max_cent) then
write(LuOut,*) ' geom_atomct_set: too many centers ',ncent,
$ names(geom)(1:lenn(geom))
geom_atomct_set = .false.
return
end if
c
do i = 1, ncent
geom_atomct(i,geom) = atomct(i)
enddo
c
end
C>
C> \brief Retrieve the atom constraint type of the centers in a geometry
C> instance
C>
C> \return Return .true. if successfull, and .false. otherwise
C>
logical function geom_atomct_get(geom, ncent, atomct)
implicit none
#include "nwc_const.fh"
#include "geomP.fh"
#include "stdio.fh"
c
integer geom !< [Input] the geometry handle
integer ncent !< [Output] the number of centers
double precision atomct(ncent) !< [Output] the atom constraint
!< type on each center
c
integer i
c
logical geom_check_handle
external geom_check_handle
c
geom_atomct_get = geom_check_handle(geom, 'geom_atomct_set')
if (.not. geom_atomct_get) return
c
ncent = ncenter(geom)
c
do i = 1, ncent
atomct(i) = geom_atomct(i,geom)
enddo
c
end
c
C> \brief Set the Angstrom to Bohr conversion factor for a geometry
C> instance

3
src/geom/geom.fh
View file

@ -72,6 +72,9 @@ c
logical geom_mass_set
logical geom_mass_get
logical geom_tag_to_covalent_radius
c
logical geom_atomct_set
logical geom_atomct_get
c
logical geom_nucexps_set
logical geom_nucexps_get

7
src/geom/geomP.fh
View file

@ -28,6 +28,9 @@ c coords(1:3,1:max_cent,1:max_geom) = cartesian coords of this geometry
c charge(1:max_cent,1:max_geom) = charges associated with centers
c geom_mass(1:max_cent,1:max_geom) = mass associated with centers
c geom_invnucexp(1:max_cent,1:max_geom) = inverse of nuclear exponent
c geom_atomct(1:max_cent,1:max_geom) = atom constraint type (keeps track internal
c coordinate constraints and their effect on
c the symmetry)
c dipole ... not yet
c quadrupole ... not yet
c pseudopotential ... not yet
@ -138,7 +141,8 @@ c
$ recip_lat_angles, sym_ops, angstrom_to_au, geom_invnucexp,
$ geom_mass,bmatrix, velocities, zmt_varsign, zmt_cvr_scaling,
$ zmt_ijbond_val, zmt_ijkang_val, zmt_ijklto_val,
$ zmt_ijklop_val, zmt_ijklnb_val, zmt_izfrz_val
$ zmt_ijklop_val, zmt_ijklnb_val, zmt_izfrz_val,
$ geom_atomct
common/cgeometry/
* doubles
@ -146,6 +150,7 @@ c
$ velocities(3,max_cent,max_geom),
$ charge(max_cent,max_geom),
$ geom_mass(max_cent,max_geom),
$ geom_atomct(max_cent,max_geom),
$ geom_invnucexp(max_cent,max_geom),
$ efield(3, max_geom),
$ latvec(3,3, max_geom), ! Obsolete ?

10
src/geom/geom_getsym.F
View file

@ -3,7 +3,7 @@ C> @{
subroutine geom_getsym0(rtdb,geom,oautosym,
, ncenter,isys,
, name,tags,
, threquiv,scale,coords,charge,velocities)
, threquiv,scale,coords,charge,velocities,atomct)
C $Id$
implicit none
#include "errquit.fh"
@ -21,7 +21,7 @@ c
integer isys
integer ncenter
double precision scale,coords(3,*),charge(*),
, velocities(*),threquiv
, velocities(*),threquiv,atomct(*)
c
logical oprint,oprint_sym
integer ncenter_before_auto
@ -43,8 +43,8 @@ c
$ call errquit('geom_input: geom_vel_get failed', 0,
& GEOM_ERR)
ncenter_before_auto = ncenter
call geom_auto_sym(rtdb,geom,coords,charge,tags,ncenter,
$ threquiv,groupname,velocities)
call geom_auto_sym(rtdb,geom,coords,charge,tags,atomct,
$ ncenter,threquiv,groupname,velocities)
if (geom_group_set(geom,groupname)) then
if (.not.geom_cart_set(geom,ncenter,tags,coords,charge))
$ call errquit('geom_input: geom_cart_set failed', 0,
@ -125,7 +125,7 @@ c
, ncenter(geom),isystype(geom),
, name,tags(1,geom),
, threquiv,scale,coords(1,1,geom),
, charge(1,geom),velocities(1,1,geom))
, charge(1,geom),velocities(1,1,geom),geom_atomct(1,geom))
c
c store stuff to rtdb
c

198
src/geom/geom_input.F
View file

@ -49,6 +49,40 @@ c . and for forcing symmetry
double precision charge(max_center), mass(max_center)
double precision invnucexp(max_center)
double precision ainv(3,3),t(3)
c
c Constraints in the geometry change the types of the atoms
c involved. E.g. in an H2O molecule the H involved in a bond of
c constraint length is not of the same type as the other one. See
c test case h2o_opt in the QA for an example. There are three
c different constraints:
c
c bond lenghts: A - A
c bond angles: B - C - B
c torsion angles: D - E - E - D
c
c These types can be represented as floating point values
c
c A: sqrt(2) B: sqrt(3) C: sqrt(5)
c D: sqrt(7) E: sqrt(11)
c
c In the code they are represented by names:
c
c A: bond_constraint_type
c B: end_angle_constraint_type
c C: center_angle_constraint_type
c D: end_torsion_constraint_type
c E: center_torsion_constraint_type
c
c This ensures that n*X .ne. m*Y for all non-zero integer n and m,
c and all types X and Y when X .ne. Y. The atom constraint types
c will be stored in the ATOMCT array. This array must be updated for
c every CONSTANT line in a zmatrix or zcoord block. The values are
c computed in the GEOM_ZMT_GEO and HND_AUTOZ_INPUT subroutines and
c used in the HND_MOLOPS subroutine in the SAMEATM statement
c function.
c
double precision atomct(max_center) !< Atom constraint type
c
character*16 tags(max_center)
character*2 symbol
character*16 element
@ -237,6 +271,9 @@ c
if (.not. geom_rtdb_load(rtdb, geom, name)) call errquit
$ ('geom_input: cannot find geometry to adjust',0,
& RTDB_ERR)
if (.not. geom_atomct_get(geom,ncenter,atomct)) call errquit
$ ('geom_input: cannot find current atom constraint types',
$ geom,GEOM_ERR)
if (oprint) write(LuOut,778) name(1:inp_strlen(name))
778 format(/,' Adjusting existing geometry named ',a,/)
endif
@ -363,6 +400,18 @@ c if(.not.geom_create_from_file(filename,rtdb)) call errquit
oautoz = .false.
goto 20
else if (inp_compare(.false.,'zcoord', field)) then
c
c When the coordinates are being changed with ZCOORD the code
c should not enforce a previously found point group. The
c ZCOORD input very likely is going to change that point
c group. So, delete the symmetry information. (Of course we
c need to re-detect the symmetry when we are done (if needed)
c but that happens automatically).
c
if (.not.geom_strip_sym(geom)) then
call errquit('geom_input: zcoord: could not delete '
$ //'symmetry information',geom,GEOM_ERR)
endif
call geom_autoz_input(geom,oprint)
goto 20
else if (inp_compare(.false.,'zmt', field).or.
@ -375,7 +424,7 @@ c if(.not.geom_create_from_file(filename,rtdb)) call errquit
if (oadjust) call errquit
$ ('geom_input: zmatrix not allowed for adjust',0,
& GEOM_ERR)
call geom_zmt_input(geom,coords,tags,charge,mass,
call geom_zmt_input(geom,coords,tags,charge,mass,atomct,
$ ncenter,scale,oprint,found_cart)
c
c If (ofinite = .true.) set invnucexp finite nucleus values
@ -572,6 +621,10 @@ c
& call errquit('geom_input: geom_nucexps_set failed', 0,
& GEOM_ERR)
if (.not. geom_atomct_set(geom, ncenter, atomct))
& call errquit('geom_input: geom_atomct_set failed', 0,
& GEOM_ERR)
c
c Check to see if both autosym and symmetry are being used. If so,
@ -590,7 +643,7 @@ c
call geom_getsym0(rtdb,geom,oautosym,
, ncenter,isys,
, name,tags,
, threquiv,scale,coords,charge,velocities)
, threquiv,scale,coords,charge,velocities,atomct)
endif ! End of else clause for if(oadjust)
c
@ -1890,6 +1943,7 @@ c
$ zmt_ijklto_frz(1,geom),
$ zmt_ijklop_frz(1,geom),
$ zmt_ijklnb_frz(1,geom),
$ geom_atomct(1,geom),
$ max_zcoord,
$ zmt_cvr_scaling(geom),
$ zmt_maxtor(geom), oprint, scale)
@ -1900,7 +1954,7 @@ c
$ ijbond_nam,ijkang_nam,ijklto_nam,ijklop_nam,ijklnb_nam,
$ ijbond_val,ijkang_val,ijklto_val,ijklop_val,ijklnb_val,
$ ijbond_frz,ijkang_frz,ijklto_frz,ijklop_frz,ijklnb_frz,
$ max_zcoord, cvr_scaling,maxtor,oprint,scale)
$ atomct, max_zcoord, cvr_scaling,maxtor,oprint,scale)
implicit none
#include "errquit.fh"
#include "inp.fh"
@ -1918,6 +1972,7 @@ c
$ ijklop_nam(*),ijklnb_nam(*)
double precision ijbond_val(*),ijkang_val(*),ijklto_val(*),
$ ijklop_val(*),ijklnb_val(*)
double precision atomct(*)
logical ijbond_frz(*),ijkang_frz(*),ijklto_frz(*),ijklop_frz(*),
$ ijklnb_frz(*)
character*13 fff
@ -1926,6 +1981,17 @@ c
logical oprint
c
integer nbond, nangl, ntors, nopla, nlinb, i, j, nuser
c
double precision bond_constraint_type
double precision end_angle_constraint_type
double precision center_angle_constraint_type
double precision end_torsion_constraint_type
double precision center_torsion_constraint_type
parameter(bond_constraint_type = sqrt(2.0d0))
parameter(end_angle_constraint_type = sqrt(3.0d0))
parameter(center_angle_constraint_type = sqrt(5.0d0))
parameter(end_torsion_constraint_type = sqrt(7.0d0))
parameter(center_torsion_constraint_type = sqrt(11.0d0))
c
cvr_scaling=cvr_factor
nbond=0
@ -2090,7 +2156,42 @@ c$$$ $ ('hnd_autoz_input: reading ijklnb',4)
endif
goto 10 ! Read next line
c
999 if (cvr_scaling.ne.cvr_factor .and. oprint)
999 continue
c
c Update the atom constraint type array
c
do i = 1, nbond
if (ijbond_frz(i)) then
atomct(ijbond(1,i)) = atomct(ijbond(1,i))
+ + bond_constraint_type
atomct(ijbond(2,i)) = atomct(ijbond(2,i))
+ + bond_constraint_type
endif
enddo
do i = 1, nangl
if (ijkang_frz(i)) then
atomct(ijkang(1,i)) = atomct(ijkang(1,i))
+ + end_angle_constraint_type
atomct(ijkang(2,i)) = atomct(ijkang(2,i))
+ + center_angle_constraint_type
atomct(ijkang(3,i)) = atomct(ijkang(3,i))
+ + end_angle_constraint_type
endif
enddo
do i = 1, ntors
if (ijklto_frz(i)) then
atomct(ijklto(1,i)) = atomct(ijklto(1,i))
+ + end_torsion_constraint_type
atomct(ijklto(2,i)) = atomct(ijklto(2,i))
+ + center_torsion_constraint_type
atomct(ijklto(3,i)) = atomct(ijklto(3,i))
+ + center_torsion_constraint_type
atomct(ijklto(4,i)) = atomct(ijklto(4,i))
+ + end_torsion_constraint_type
endif
enddo
c
if (cvr_scaling.ne.cvr_factor .and. oprint)
$ write(LuOut,9) cvr_scaling
9 format(' cvr_scaling: ', f6.2)
if (maxtor.ne.100 .and. oprint) write(LuOut,5) maxtor
@ -2140,7 +2241,7 @@ c
c
end
subroutine geom_zmt_input(geom,
1 coords,tags,charge,mass,ncenter,
1 coords,tags,charge,mass,atomct,ncenter,
2 units,oprint,found_cart)
implicit none
#include "errquit.fh"
@ -2197,6 +2298,7 @@ c
double precision coords(3,*)
double precision charge(*)
double precision mass(*)
double precision atomct(*)
dimension tags(*)
c
iw = luout
@ -2210,7 +2312,7 @@ c
c
c ----- call routine to read z-matrix data -----
c
call geom_zmt_geo(coords,charge,tags,ncenter,
call geom_zmt_geo(coords,charge,tags,atomct,ncenter,
1 nizmat,izmat,nzvar,izfrz,nzfrz,
2 units,zvarname,zvarsign,found_cart)
c
@ -2920,7 +3022,7 @@ c
enddo
c
end
SUBROUTINE GEOM_ZMT_GEO(COORDS,CHARGE,TAGS,NCENTER,
SUBROUTINE GEOM_ZMT_GEO(COORDS,CHARGE,TAGS,ATOMCT,NCENTER,
1 IZ,IZMAT,NZMOD,ICFRZ,NCFRZ,
2 UNITS,zvarname,zvarsign,
3 found_cart)
@ -3029,6 +3131,7 @@ C
1 ATNAME(MAXGEO),ATNUM(MAXGEO),CART(MAXGEO)
double precision COORDS(3,*)
double precision CHARGE( *)
double precision ATOMCT( *)
DIMENSION TAGS( *)
DIMENSION IZMAT( *)
DIMENSION ICFRZ( *)
@ -3064,6 +3167,16 @@ C
logical inp_compare
double precision zero
double precision xxiat,yyiat,zziat
double precision bond_constraint_type
double precision end_angle_constraint_type
double precision center_angle_constraint_type
double precision end_torsion_constraint_type
double precision center_torsion_constraint_type
parameter(bond_constraint_type = sqrt(2.0d0))
parameter(end_angle_constraint_type = sqrt(3.0d0))
parameter(center_angle_constraint_type = sqrt(5.0d0))
parameter(end_torsion_constraint_type = sqrt(7.0d0))
parameter(center_torsion_constraint_type = sqrt(11.0d0))
EQUIVALENCE (BLNK80,BLK(1))
EQUIVALENCE (CHREND,WRDEND)
DATA ERRMSG /'PROGRAM ','STOP IN ','- ZGEO -'/
@ -3287,6 +3400,7 @@ C ATOMS IS DETECTED VIA A BLANK LINE
C
IAT = 0
110 IAT = IAT + 1
ATOMCT(IAT) = 0.0d0
CALL HND_RDFREE(IR,STRING,IERR)
IF(IERR.NE.0) THEN
WRITE(IW,8886)
@ -3471,6 +3585,51 @@ C
C ----- WE HAVE THE -Z- MATRIX WITH THE INDICES IN -ZMT- -----
C AND THE VALUES IN -ZVAL-
C
c -FRZVAL- lists which constraints are in place
c FRZVAL(1,iat).eqv..true. : bond length is a constant
c FRZVAL(2,iat).eqv..true. : bond angle is a constant
c FRZVAL(3,iat).eqv..true. : torsion angle is a constant
c Combined with the indices in -ZMT- this is all we need to
c compute -ATOMCT-.
c
c update -ATOMCT- for constant bond lengths
c
do iat = 2, nat
if (frzval(1,iat)) then
atomct(zmt(1,iat)) = atomct(zmt(1,iat))
+ + bond_constraint_type
atomct(zmt(2,iat)) = atomct(zmt(2,iat))
+ + bond_constraint_type
endif
enddo
c
c update -ATOMCT- for constant bond angles
c
do iat = 3, nat
if (frzval(2,iat)) then
atomct(zmt(1,iat)) = atomct(zmt(1,iat))
+ + end_angle_constraint_type
atomct(zmt(2,iat)) = atomct(zmt(2,iat))
+ + center_angle_constraint_type
atomct(zmt(3,iat)) = atomct(zmt(3,iat))
+ + end_angle_constraint_type
endif
enddo
c
c update -ATOMCT- for constant torsion angles
c
do iat = 4, nat
if (frzval(3,iat)) then
atomct(zmt(1,iat)) = atomct(zmt(1,iat))
+ + end_torsion_constraint_type
atomct(zmt(2,iat)) = atomct(zmt(2,iat))
+ + center_torsion_constraint_type
atomct(zmt(3,iat)) = atomct(zmt(3,iat))
+ + center_torsion_constraint_type
atomct(zmt(4,iat)) = atomct(zmt(4,iat))
+ + end_torsion_constraint_type
endif
enddo
IF(DBUG) THEN
WRITE(IW,9998) TITLE
WRITE(IW,9997)
@ -3482,6 +3641,7 @@ C
WRITE(IW,9996) IAT,( ZMT(J,IAT),J=1,4),
1 ( ZVAL(J,IAT),J=1,3),
2 (FRZVAL(J,IAT),J=1,3)
write(iw,9970) IAT, atomct(iat)
ENDIF
ENDDO
IF(LST) THEN
@ -3746,6 +3906,7 @@ C
9976 FORMAT(' ----- -LST- POINT NO. = ',I3,' -----')
9975 FORMAT(' $PES NPES =',I4,', IUNIT =1, $END')
9974 FORMAT(' IAT=',I5,' ATNAME=',2X,A8)
9970 FORMAT(i5,' ATOMCT=',f15.5)
8889 FORMAT(/,10X,12(1H-),/,10X,'-$GEO- INPUT',/,10X,12(1H-))
8888 FORMAT(A8)
8887 FORMAT(' NO DATA GROUP -$GEO- FOUND AS ALTERNATE INPUT. ')
@ -3986,8 +4147,9 @@ c
& 1x,' which is larger than the threshold:',6x,f20.13,/,/)
end
subroutine geom_auto_sym(rtdb,geom,
1 coords,charge,tags,ncenter,threquiv,
$ group,veloct)
1 coords,charge,tags,atomct,
& ncenter,threquiv,
$ group,veloct)
implicit none
#include "stdio.fh"
#include "nwc_const.fh"
@ -4010,6 +4172,7 @@ c
double precision v
double precision veloct
double precision zan
double precision atomct(*)
character*16 tags
character*16 atmlab
character*8 groupname
@ -4063,7 +4226,7 @@ c
c
c ----- call auto_sym -----
c
call hnd_autsym(odone,rtdb,threquiv,groupname,geom)
call hnd_autsym(atomct,odone,rtdb,threquiv,groupname,geom)
c
group = groupname
c
@ -4082,7 +4245,7 @@ c
9997 format(1x,i5,i5,3f10.6,f10.3)
c
end
SUBROUTINE HND_AUTSYM(ODONE,RTDB,THREQUIV,groupname,geom)
SUBROUTINE HND_AUTSYM(atomct,ODONE,RTDB,THREQUIV,groupname,geom)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
implicit none
#include "nwc_const.fh"
@ -4113,6 +4276,7 @@ c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
double precision AXS(3,3),EIG(3)
double precision RT(3,3)
double precision TR(3)
double precision atomct(*)
integer iat,i
C
DBUG=.FALSE.
@ -4154,8 +4318,8 @@ C
C
C ----- IDENTIFY SYMMETRY OPERATIONS -----
C
CALL HND_MOLOPS(C,C1,C2,C3,NAT,AXS,EIG,RT,TR,ODONE,THREQUIV,
$ groupname,geom,V,V2,V3)
CALL HND_MOLOPS(C,C1,C2,C3,atomct,NAT,AXS,EIG,RT,TR,ODONE,
$ THREQUIV,groupname,geom,V,V2,V3)
C
C ----- CREATE ATOM AND SHELL MAPPINGS -----
C
@ -4192,7 +4356,7 @@ C
2 /,10X,7(1H-))
9996 FORMAT(/)
END
SUBROUTINE HND_MOLOPS(C,C1,C2,C3,NAT,AXS,EIG,RT,TR,ODONE,
SUBROUTINE HND_MOLOPS(C,C1,C2,C3,atomct,NAT,AXS,EIG,RT,TR,ODONE,
$ THREQUIV,groupname,geom,V,V2,V3)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
implicit none
@ -4203,6 +4367,8 @@ C CODE ADAPTED FROM A.M.CHAKA's ORIGINAL CODE.
c
c RJH ... added parameter THREQUIV to provide user control
c . for noisy geometries, and groupname to return value
c HvD ... Added ATOMCT to allow checking for symmetry breaking
c constraints on the geometry
C
#include "util.fh"
#include "nwc_const.fh"
@ -4281,6 +4447,7 @@ C
double precision V3(3,*)
double precision CM(3),AM(3,3)
double precision AXM(3,3)
double precision atomct(*)
logical C2AXS(3)
logical C4AXS(3)
logical S4AXS(3)
@ -4314,7 +4481,8 @@ c DATA TENM02 /1.0D-02/
DATA TENM04 /1.0D-04/
DATA TENM05 /1.0D-05/
DATA EPS /1.01D+00/
sameatm(JAT,IAT)=ATMLAB(JAT).EQ.ATMLAB(IAT)
sameatm(JAT,IAT)=ATMLAB(JAT).EQ.ATMLAB(IAT).and.
+ (abs(atomct(jat)-atomct(iat)).lt.1.0d-8)
big=dlamch('o')**(0.125d0)
QZERO=cmplx(ZERO,ZERO)
PI= FOUR*ATAN(ONE)