Line data Source code
1 : !!****m* ABINIT/m_abimover
2 : !! NAME
3 : !! m_abimover
4 : !!
5 : !! FUNCTION
6 : !! This module contains definition the types abimover, mttk, abiforstr, delocint, and bonds
7 : !! and their related ini and free routines
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2001-2026 ABINIT group (DCA, XG, GMR, SE, Mver, JJ)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_abimover
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_atomdata
28 : use m_errors
29 : use m_dtset
30 : use m_dtfil
31 :
32 : use m_geometry, only : acrossb
33 : !use m_fstrings, only : sjoin, itoa
34 :
35 : implicit none
36 :
37 : private
38 :
39 : public :: abimover_ini
40 : public :: abimover_destroy
41 : public :: mttk_ini ! initialize the object
42 : public :: mttk_fin ! Release memory
43 : public :: abiforstr_ini ! Initialize the object
44 : public :: abiforstr_fin ! Free memory
45 : public :: delocint_ini ! Initialize the delocint object
46 : public :: delocint_fin ! Free memory
47 : public :: bonds_free
48 : public :: bond_length
49 : public :: print_bonds
50 : public :: make_bonds_new
51 : public :: calc_prim_int
52 : public :: make_prim_internals
53 :
54 : integer,public, parameter :: mover_BEFORE=0
55 : integer,public, parameter :: mover_AFTER=1
56 : !!***
57 :
58 : !----------------------------------------------------------------------
59 :
60 : !!****t* m_abimover/abimover
61 : !! NAME
62 : !! abimover
63 : !!
64 : !! FUNCTION
65 : !! This datatype has the purpose of store all the data taken
66 : !! usually from dtset (but not only) needed for the different predictors
67 : !! to update positions, acell, etc.
68 : !!
69 : !! SOURCE
70 :
71 : type, public :: abimover
72 :
73 : ! scalars
74 : ! Delay of Permutation (Used by pred_langevin only)
75 : integer :: delayperm
76 : ! DIIS memory (Used by pred_diisrelax only)
77 : integer :: diismemory
78 : ! Geometry Optimization Precondition option
79 : integer :: goprecon
80 : ! include a JELLium SLAB in the cell
81 : integer :: jellslab
82 : ! Number of ATOMs
83 : integer :: natom
84 : ! Number of CONstraint EQuations
85 : integer :: nconeq
86 : ! Number of Degrees Of Freedom
87 : integer :: ndof
88 : ! number of Shifts for the Qpoint Grid (used for ionmov 26 and 27)
89 : integer :: ph_nqshift
90 : ! Use by pred_isothermal only
91 : integer :: nnos
92 : ! Number of SYMmetry operations
93 : integer :: nsym
94 : ! Number of Types of atoms
95 : integer :: ntypat
96 : ! OPTimize the CELL shape and dimensions
97 : integer :: optcell
98 : ! RESTART Xcart and Fred
99 : integer :: restartxf
100 : ! Sign of Permutation (Used by pred_langevin only)
101 : integer :: signperm
102 : ! Ion movement
103 : integer :: ionmov
104 :
105 : ! Flag: true if we use a PIMD version of the mover routine
106 : logical :: use_pimd_routine
107 :
108 : ! Use by pred_isothermal only
109 : real(dp) :: bmass
110 : ! Delta Time for IONs
111 : real(dp) :: dtion
112 : ! Used by pred_langevin only
113 : real(dp) :: friction
114 : ! Used by pred_langevin only
115 : real(dp) :: mdwall
116 : ! Used by pred_nose only
117 : real(dp) :: noseinert
118 : ! STRess PRECONditioner
119 : real(dp) :: strprecon
120 : ! VIScosity
121 : real(dp) :: vis
122 :
123 : ! arrays
124 : ! Indices of AToms that are FIXed
125 : integer,pointer :: iatfix(:,:) ! iatfix(3,natom)
126 : ! SYMmetries, Anti-FerroMagnetic characteristics
127 : integer,pointer :: symafm(:) ! symafm(nsym)
128 : ! SYMmetry in REaL space
129 : integer,pointer :: symrel(:,:,:) ! symrel(3,3,nsym)
130 : ! Translation NON-Symmorphic vectors
131 : real(dp),pointer :: tnons(:,:) ! tnons(3,nsym)
132 : ! TYPe of ATom
133 : integer,pointer :: typat(:) ! typat(natom)
134 : ! PRTint ATom LIST
135 : integer,pointer :: prtatlist(:) ! prtatlist(natom)
136 : ! Qpoint grid (used for ionmov 26 and 27)
137 : integer,pointer :: ph_ngqpt(:) ! ph_ngqpt(3)
138 : ! shift of the Qpoint Grid (used for ionmov 26 and 27)
139 : real(dp),pointer :: ph_qshift(:,:) !
140 : ! amu input var for the current image
141 : real(dp), pointer :: amu_curr(:) ! amu_curr(ntypat)
142 : ! Mass of each atom
143 : real(dp),pointer :: amass(:) ! amass(natom)
144 : ! Geometry Optimization Preconditioner PaRaMeters
145 : real(dp),pointer :: goprecprm(:)
146 : ! Molecular Dynamics Initial and Final Temperature
147 : real(dp),pointer :: mdtemp(:) ! mdtemp(2) (initial,final)
148 : ! STRess TARGET
149 : real(dp),pointer :: strtarget(:) ! strtarget(6)
150 : ! Use by pred_isothermal only
151 : real(dp),pointer :: qmass(:)
152 : ! Z number of each NUCLeus
153 : real(dp),pointer :: znucl(:) ! znucl(npsp)
154 :
155 : ! Filename for Hessian matrix
156 : character(len=fnlen), pointer :: fnameabi_hes
157 : ! Filename for _HIST file
158 : character(len=fnlen), pointer :: filnam_ds(:) ! dtfil%filnam_ds(5)
159 :
160 : end type abimover
161 : !!***
162 :
163 : !----------------------------------------------------------------------
164 :
165 : !!****t* m_abimover/abimover_specs
166 : type,public :: abimover_specs
167 :
168 : !scalars
169 : integer :: ncycle
170 : integer :: nhist ! Number of step of history needed in the algorithm
171 : character(len=8) :: crit4xml
172 : character(len=10) :: type4xml
173 : character(len=60) :: method
174 : logical :: isFconv ! If the convergence is needed
175 : logical :: isARused
176 : logical :: isVused
177 :
178 : end type abimover_specs
179 : !!***
180 :
181 : !----------------------------------------------------------------------
182 :
183 : !!****t* m_abimover/delocint
184 : !! NAME
185 : !! delocint
186 : !!
187 : !! FUNCTION
188 : !! Datatype with the important variables in pred_delocint
189 : !!
190 : !! NOTES
191 : !! deloc <type(delocint)>=Important variables for pred_delocint
192 : !! |
193 : !! ! icenter = Index of the center of the number of shifts
194 : !! | nang = Number of angles
195 : !! | nbond = Number of bonds
196 : !! | ncart = Number of cartesian directions (used for constraints)
197 : !! | ndihed = Number of dihedrals
198 : !! | nrshift = Dimension of rshift
199 : !! | ninternal= Number of internal coordinates
200 : !! | ninternal=nbond+nang+ndihed+ncart
201 : !! | angs(2,3,nang) = Indexes to characterize angles
202 : !! | bonds(2,2,nbond)= For a bond between iatom and jatom
203 : !! | bonds(1,1,nbond) = iatom
204 : !! | bonds(2,1,nbond) = icenter
205 : !! | bonds(1,2,nbond) = jatom
206 : !! | bonds(2,2,nbond) = irshift
207 : !! | carts(2,ncart) = Index of total primitive internal, and atom (carts(2,:))
208 : !! | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals
209 : !! | rshift(3,nrshift)= Shift in xred that must be done to find
210 : !! | all neighbors of a given atom within a
211 : !! | given number of neighboring shells
212 : !!
213 : !! SOURCE
214 :
215 : type,public :: delocint
216 :
217 : ! scalars
218 : integer :: icenter
219 : integer :: nang
220 : integer :: nbond
221 : integer :: ncart
222 : integer :: ndihed
223 : integer :: nrshift
224 : integer :: ninternal
225 :
226 : ! arrays
227 : integer,allocatable :: angs(:,:,:)
228 : integer,allocatable :: bonds(:,:,:)
229 : integer,allocatable :: carts(:,:)
230 : integer,allocatable :: dihedrals(:,:,:)
231 : real(dp),allocatable :: rshift(:,:)
232 :
233 : end type delocint
234 : !!***
235 :
236 : !----------------------------------------------------------------------
237 :
238 : !!****t* m_abimover/mttk_type
239 : !! NAME
240 : !! mttk_type
241 : !!
242 : !! FUNCTION
243 : !! For Martyna et al. (TTK) reversible MD integration scheme and related data
244 : !!
245 : !! SOURCE
246 :
247 : type, public :: mttk_type
248 :
249 : real(dp) :: glogv
250 : !Logarithm of the volume
251 :
252 : real(dp) :: vlogv
253 : !Derivative of logv
254 :
255 : real(dp) :: gboxg(3,3)
256 : !Imbalance in pressure (see paper)
257 :
258 : real(dp) :: vboxg(3,3)
259 : !Velocity of log rprimd (see paper)
260 :
261 : real(dp), allocatable :: glogs(:)
262 : ! glogs(nnos)
263 : ! Imbalance of kinetic energy
264 :
265 : real(dp), allocatable :: vlogs(:)
266 : ! vlogs(nnos)
267 : ! Velocities of thermostat variables
268 :
269 : real(dp), allocatable :: xlogs(:)
270 : ! xlogs(nnos)
271 : ! Positions of thermostat variables
272 :
273 : end type mttk_type
274 : !!***
275 :
276 : !----------------------------------------------------------------------
277 :
278 : !!****t* m_abimover/abiforstr
279 : !! NAME
280 : !! abiforstr
281 : !!
282 : !! FUNCTION
283 : !! Store forces, stress and energy, cartesian and reduced forces
284 : !! one scalar for energy and 6 element array for stress
285 : !!
286 : !! NOTES
287 : !!
288 : !! SOURCE
289 :
290 : type, public :: abiforstr
291 :
292 : ! scalars
293 : real(dp) :: etotal
294 : ! Total energy
295 :
296 : ! arrays
297 : real(dp),allocatable :: fcart(:,:)
298 : ! Cartesian forces
299 : real(dp),allocatable :: gred(:,:)
300 : ! Reduced forces
301 : real(dp) :: strten(6)
302 : ! Stress tensor (Symmetrical 3x3 matrix)
303 :
304 : end type abiforstr
305 : !!***
306 :
307 : !----------------------------------------------------------------------
308 :
309 : !!****t* m_abimover/ab_xfh_type
310 : !! NAME
311 : !! ab_xfh_type
312 : !!
313 : !! FUNCTION
314 : !! Datatype with the old structure for storing history
315 : !! used in gstate and brdmin,delocint, and others
316 : !!
317 : !! NOTES
318 : !! This is a transitional structure, to bridge between
319 : !! the old code and the new one base on abihist
320 : !!
321 : !! SOURCE
322 :
323 : type, public :: ab_xfh_type
324 :
325 : integer :: nxfh,nxfhr,mxfh
326 : ! mxfh = last dimension of the xfhist array
327 : ! nxfh = actual number of (x,f) history pairs, see xfhist array
328 :
329 : real(dp),allocatable :: xfhist(:,:,:,:)
330 : ! xfhist(3,natom+4,2,mxfh) = (x,f) history array, also including rprim and stress
331 :
332 : end type ab_xfh_type
333 : !!***
334 :
335 : !----------------------------------------------------------------------
336 :
337 : !!****t* m_abimover/go_bonds
338 : !! NAME
339 : !! go_bonds
340 : !!
341 : !! FUNCTION
342 : !! Datatype all the information relevant to create
343 : !! bonds between atoms inside and outside the cell
344 : !!
345 : !! SOURCE
346 :
347 : type, public :: go_bonds
348 :
349 : !scalar
350 : real(dp) :: tolerance ! To decide if consider bond the atom or not
351 : ! 1.0 means that only consider values lower
352 : ! than the sum of covalent radius
353 :
354 : integer :: nbonds ! Total number of bonds for the system
355 :
356 : !arrays
357 :
358 : integer,allocatable :: nbondi(:) ! Number of bonds for atom i
359 : integer,allocatable :: indexi(:,:) ! Indices of bonds for atom i
360 : ! Positive: Vector from i to j
361 : ! Negative: Vector from j to i
362 :
363 : real(dp),allocatable :: bond_length(:) ! Bond lengths
364 : real(dp),allocatable :: bond_vect(:,:) ! Unitary vectors for bonds
365 :
366 : end type go_bonds
367 : !!***
368 :
369 : !----------------------------------------------------------------------
370 :
371 : !!****t* m_abimover/go_angles
372 : !! NAME
373 : !! go_angles
374 : !!
375 : !! FUNCTION
376 : !! Datatype all the information relevant to create
377 : !! angles between atoms inside and outside the cell
378 : !!
379 : !! NOTES
380 : !! This type is not used
381 : !!
382 : !! SOURCE
383 :
384 : type, public :: go_angles
385 :
386 : !scalar
387 : integer :: nangles ! Total number of bonds for the system
388 :
389 : !arrays
390 : integer,allocatable :: angle_vertex(:) ! Indices of the vertex atom
391 : real(dp),allocatable :: angle_value(:) ! Value of angle in radians
392 : real(dp),allocatable :: angle_bonds(:,:) ! Indices of the bonds
393 : real(dp),allocatable :: angle_vect(:,:) ! Unitary vector perpendicular to the plane
394 :
395 : end type go_angles
396 :
397 : !public :: make_angles_new ! This routine is broken and should be tested before use.
398 : !!***
399 :
400 : !----------------------------------------------------------------------
401 :
402 : contains !=============================================================
403 : !!***
404 :
405 :
406 : !!****f* m_abimover/abimover_ini
407 : !! NAME
408 : !! abimover_ini
409 : !!
410 : !! FUNCTION
411 : !! Initializes the abimover structure and the abimover_specs information
412 : !!
413 : !! INPUTS
414 : !!
415 : !! OUTPUT
416 : !!
417 : !! NOTES
418 : !!
419 : !! SOURCE
420 :
421 4368 : subroutine abimover_ini(ab_mover,amu_curr,dtfil,dtset,specs)
422 :
423 : !Arguments ------------------------------------
424 : real(dp),target, intent(in) :: amu_curr(:) ! amu_curr(ntype)
425 : type(abimover),intent(out) :: ab_mover
426 : type(datafiles_type),target,intent(in) :: dtfil
427 : type(dataset_type),target,intent(in) :: dtset
428 : type(abimover_specs),intent(out) :: specs
429 :
430 : !Local variables-------------------------------
431 : !scalars
432 : integer :: natom
433 : character(len=500) :: msg
434 : !arrays
435 :
436 : ! ***************************************************************
437 :
438 :
439 : !write(std_out,*) 'mover 01'
440 : !###########################################################
441 : !### 01. Initialization of ab_mover
442 :
443 : !Copy or create pointers for the information from the Dataset (dtset) to the ab_mover structure
444 4368 : natom=dtset%natom
445 :
446 4368 : ab_mover%delayperm =dtset%delayperm
447 4368 : ab_mover%diismemory =dtset%diismemory
448 4368 : ab_mover%goprecon =dtset%goprecon
449 4368 : ab_mover%jellslab =dtset%jellslab
450 4368 : ab_mover%natom =dtset%natom
451 4368 : ab_mover%nconeq =dtset%nconeq
452 4368 : ab_mover%nnos =dtset%nnos
453 4368 : ab_mover%nsym =dtset%nsym
454 4368 : ab_mover%ntypat =dtset%ntypat
455 4368 : ab_mover%optcell =dtset%optcell
456 4368 : ab_mover%restartxf =dtset%restartxf
457 4368 : ab_mover%signperm =dtset%signperm
458 4368 : ab_mover%ionmov =dtset%ionmov
459 4368 : ab_mover%bmass =dtset%bmass
460 4368 : ab_mover%dtion =dtset%dtion
461 4368 : ab_mover%friction =dtset%friction
462 4368 : ab_mover%mdwall =dtset%mdwall
463 4368 : ab_mover%noseinert =dtset%noseinert
464 4368 : ab_mover%ph_nqshift =dtset%ph_nqshift
465 4368 : ab_mover%strprecon =dtset%strprecon
466 4368 : ab_mover%vis =dtset%vis
467 :
468 4368 : ab_mover%iatfix =>dtset%iatfix(:,1:natom)
469 4368 : ab_mover%symafm =>dtset%symafm
470 4368 : ab_mover%symrel =>dtset%symrel
471 4368 : ab_mover%tnons =>dtset%tnons
472 4368 : ab_mover%ph_ngqpt =>dtset%ph_ngqpt
473 4368 : ab_mover%ph_qshift =>dtset%ph_qshift
474 4368 : ab_mover%typat =>dtset%typat(1:natom)
475 4368 : ab_mover%prtatlist =>dtset%prtatlist(1:natom)
476 4368 : ab_mover%goprecprm =>dtset%goprecprm
477 4368 : ab_mover%mdtemp =>dtset%mdtemp
478 4368 : ab_mover%strtarget =>dtset%strtarget
479 4368 : ab_mover%qmass =>dtset%qmass
480 4368 : ab_mover%znucl =>dtset%znucl
481 :
482 4368 : ab_mover%use_pimd_routine = .false.
483 :
484 4368 : ab_mover%amu_curr =>amu_curr
485 13104 : ABI_MALLOC(ab_mover%amass,(natom))
486 39310 : ab_mover%amass(1:natom) = amu_emass * amu_curr(dtset%typat(1:natom))
487 :
488 : !Count the number of degrees of freedom, taking into account iatfix.
489 : !Ndof = 3N (when no atom is fixed). This value may be adjusted for
490 : !each algorithm below.
491 65516 : ab_mover%ndof=count(ab_mover%iatfix==0)
492 :
493 : !Filename for Hessian matrix (NOT IN DTSET)
494 4368 : ab_mover%fnameabi_hes =>dtfil%fnameabi_hes
495 : !Filename for _HIST file
496 4368 : ab_mover%filnam_ds =>dtfil%filnam_ds
497 :
498 : !call abimover_print(ab_mover,ab_out)
499 :
500 : !write(std_out,*) 'mover 02'
501 : !###########################################################
502 : !### 02. Particularities of each predictor
503 :
504 : !Default values first
505 : !--------------------
506 :
507 : !acell and rprimd are never changed except if optcell/=0
508 4368 : if (ab_mover%optcell/=0)then
509 215 : specs%isARused=.TRUE.
510 : else
511 4153 : specs%isARused=.FALSE.
512 : end if
513 :
514 : !Velocities are never changed excepts for ionmov=1,6,7,8
515 4368 : specs%isVused=.FALSE.
516 :
517 : !In general convergence is needed
518 4368 : specs%isFconv=.TRUE.
519 :
520 : !specs%ncycle is 1 by default except for ionmov=1,9,14
521 4368 : specs%ncycle=1
522 :
523 : !specs%nhist is -1 by default store all the history
524 4368 : specs%nhist=-1
525 :
526 : !This is the initialization for ionmov==1
527 : !-----------------------------------------
528 12 : select case (ab_mover%ionmov)
529 : case (1)
530 12 : specs%ncycle=4 ! Number of internal cycles for first itime
531 12 : specs%isFconv=.FALSE. ! Convergence is not used for MD
532 12 : specs%isVused=.TRUE. ! Velocities are used
533 : ! TEMPORARLY optcell is not allow
534 12 : specs%isARused=.FALSE.
535 : ! Values use in XML Output
536 12 : specs%type4xml='moldyn'
537 12 : specs%crit4xml='none'
538 : ! Name of specs%method
539 12 : if (abs(ab_mover%vis)<=1.d-8) then
540 0 : specs%method = 'Molecular dynamics without viscosity (vis=0)'
541 : else
542 : write(specs%method,'(a,1p,e12.5,a)')&
543 12 : 'Molecular dynamics with viscosity (vis=',ab_mover%vis,')'
544 : end if
545 : ! Number of history
546 12 : specs%nhist = 6
547 : ! This is the initialization for ionmov==2,3
548 : ! -------------------------------------------
549 : case (2,3)
550 : ! Values use in XML Output
551 576 : specs%type4xml='bfgs'
552 576 : specs%crit4xml='tolmxf'
553 : ! Name of specs%method
554 576 : if (ab_mover%ionmov==2) then
555 323 : specs%method = 'Broyden-Fletcher-Goldfarb-Shanno method (forces)'
556 : else
557 253 : specs%method = 'Broyden-Fletcher-Goldfarb-Shanno method (forces,Tot energy)'
558 : end if
559 : ! Number of history
560 576 : specs%nhist = 3
561 : ! This is the initialization for ionmov==4,5
562 : ! -------------------------------------------
563 : case (4,5)
564 : ! Values used in XML Output
565 14 : specs%type4xml='simple'
566 14 : specs%crit4xml='tolmxf'
567 : ! Name of specs%method
568 14 : if (ab_mover%ionmov==4) then
569 12 : specs%method = 'Conjugate gradient of potential and ionic degrees of freedom'
570 : else
571 2 : specs%method = 'Simple relaxation of ionic positions'
572 : end if
573 : ! Number of history
574 14 : specs%nhist = 3
575 : ! This is the initialization for ionmov==6
576 : ! ------------------------------------------
577 : case (6)
578 4 : specs%isFconv=.FALSE. ! Convergence is not used for MD
579 : ! TEMPORARLY optcell is not allow
580 4 : specs%isARused=.FALSE.
581 4 : specs%isVused=.TRUE. ! Velocities are used
582 : ! Values use in XML Output
583 4 : specs%type4xml='verlet'
584 4 : specs%crit4xml='tolmxf'
585 : ! Name of specs%method
586 4 : specs%method = 'Verlet algorithm for molecular dynamics'
587 : ! Number of history
588 4 : specs%nhist = 3
589 : ! This is the initialization for ionmov==7
590 : ! ------------------------------------------
591 : case (7)
592 : ! TEMPORARLY optcell is not allow
593 6 : specs%isARused=.FALSE.
594 6 : specs%isVused=.TRUE. ! Velocities are used
595 : ! Values use in XML Output
596 6 : specs%type4xml='verlet'
597 6 : specs%crit4xml='tolmxf'
598 : ! Name of specs%method
599 6 : specs%method = 'Verlet algorithm blocking every atom where dot(vel,force)<0'
600 : ! Number of history
601 6 : specs%nhist = 3
602 : ! This is the initialization for ionmov==8
603 : ! ------------------------------------------
604 : case (8)
605 4 : specs%isVused=.TRUE.
606 : ! TEMPORARLY optcell is not allow
607 4 : specs%isARused=.FALSE.
608 : ! Values use in XML Output
609 4 : specs%type4xml='nose'
610 4 : specs%crit4xml='tolmxf'
611 : ! Name of specs%method
612 4 : specs%method = 'Verlet algorithm with a nose-hoover thermostat'
613 : ! Number of history
614 4 : specs%nhist = 3
615 : ! This is the initialization for ionmov==9
616 : ! ------------------------------------------
617 : case (9)
618 : ! TEMPORARLY optcell is not allow
619 4 : specs%isARused=.FALSE.
620 4 : specs%isVused=.TRUE. ! Velocities are used
621 4 : specs%ncycle=3
622 : ! Values use in XML Output
623 4 : specs%type4xml='langevin'
624 4 : specs%crit4xml='tolmxf'
625 : ! Name of specs%method
626 4 : specs%method = 'Langevin molecular dynamics'
627 : ! Number of history
628 4 : specs%nhist = 3
629 : ! This is the initialization for ionmov==10 and 11
630 : ! -------------------------------------------
631 : case (10,11)
632 : ! TEMPORARLY optcell is not allow
633 2 : specs%isARused=.FALSE.
634 : ! Values use in XML Output
635 2 : if(ab_mover%ionmov==10)specs%type4xml='delocint'
636 2 : if(ab_mover%ionmov==11)specs%type4xml='cg'
637 2 : specs%crit4xml='tolmxf'
638 : ! Name of specs%method
639 2 : if(ab_mover%ionmov==10)specs%method = 'BFGS with delocalized internal coordinates'
640 2 : if(ab_mover%ionmov==11)specs%method = 'Conjugate gradient with deloc. int. coord.'
641 : ! Number of history
642 2 : specs%nhist = 3
643 : ! This is the initialization for ionmov==12
644 : ! -------------------------------------------
645 : case (12)
646 : ! TEMPORARLY optcell is not allow
647 50 : specs%isARused=.FALSE.
648 50 : specs%isVused=.TRUE. ! Velocities are used
649 : ! Values use in XML Output
650 50 : specs%isFconv=.FALSE. ! Convergence is not used for MD
651 50 : specs%type4xml='isokin'
652 50 : specs%crit4xml='tolmxf'
653 : ! Name of specs%method
654 50 : specs%method = 'Isokinetic ensemble molecular dynamics'
655 : ! Number of history
656 50 : specs%nhist = 3
657 : ! Number of degrees of freedom: 3N-4
658 50 : ab_mover%ndof=ab_mover%ndof-1 ! Kinetic energy conservation
659 : ! Conservation of the total momentum for each dimension, in case
660 : ! no atom position is fixed for that dimension
661 3248 : ab_mover%ndof=ab_mover%ndof-count(sum(ab_mover%iatfix,dim=2)==0)
662 : ! This is the initialization for ionmov==13
663 : ! -------------------------------------------
664 : case (13)
665 : ! optcell is allow
666 34 : specs%isARused=.TRUE. ! RPRIMD and ACELL may change
667 34 : specs%isVused=.TRUE. ! Velocities are used
668 34 : specs%isFconv=.FALSE. ! Convergence is not used for MD
669 : ! Values use in XML Output
670 34 : specs%type4xml='isother'
671 34 : specs%crit4xml='tolmxf'
672 : ! Name of specs%method
673 34 : specs%method = 'Isothermal/isenthalpic ensemble molecular dynamics'
674 : ! Number of history
675 34 : specs%nhist = 3
676 : ! This is the initialization for ionmov==14
677 : ! -------------------------------------------
678 : case (14)
679 2 : specs%ncycle=16
680 2 : specs%isFconv=.FALSE. ! Convergence is not used for MD
681 2 : specs%isVused=.TRUE. ! Velocities are used
682 : ! TEMPORARLY optcell is not allow
683 2 : specs%isARused=.FALSE.
684 : ! Values use in XML Output
685 2 : specs%type4xml='srkna14'
686 2 : specs%crit4xml='tolmxf'
687 : ! Name of specs%method
688 2 : specs%method = 'Symplectic algorithm Runge-Kutta-Nystrom SRKNa14'
689 : ! Number of history
690 2 : specs%nhist = 3
691 :
692 : ! This is the initialization for ionmov==15
693 : ! -------------------------------------------
694 : case (15)
695 : ! Values use in XML Output
696 2 : specs%type4xml='FIRE'
697 2 : specs%isVused=.TRUE. ! Velocities are used
698 2 : specs%isARused=.TRUE.
699 2 : specs%crit4xml='tolmxf'
700 : ! Name of specs%method
701 2 : specs%method = 'Fast inertial relaxation engine'
702 : ! Number of history
703 2 : specs%nhist = 2
704 : ! This is the initialization for ionmov==16
705 : ! ------------------------------------------
706 : case (16)
707 4 : ab_mover%use_pimd_routine = .true.
708 : ! TEMPORARLY optcell is not allowed
709 4 : specs%isVused=.TRUE. ! Velocities are used
710 4 : specs%isFconv=.FALSE. ! Convergence is not used for MD
711 : specs%ncycle=1
712 : ! Values use in XML Output
713 4 : specs%type4xml='langevin'
714 4 : specs%crit4xml='tolmxf'
715 : ! Name of specs%method
716 4 : specs%method = 'Langevin molecular dynamics'
717 : ! Number of history
718 4 : specs%nhist = 3
719 : ! Number of degrees of freedom: 3N-constraints
720 : ! For now the pimd langevin algorithms do not support iatfix.
721 : if(dtset%pitransform==1.or.dtset%pitransform==2.or.&
722 4 : & dtset%pimd_constraint==1.or.dtset%optcell==2) then
723 2 : ab_mover%ndof=ab_mover%ndof-3
724 : end if
725 : ! This is the initialization for ionmov==20
726 : ! -------------------------------------------
727 : case (20)
728 : ! TEMPORARLY optcell is not allow
729 0 : specs%isARused=.FALSE.
730 : ! Values use in XML Output
731 0 : specs%type4xml='diisrelax'
732 0 : specs%crit4xml='tolmxf'
733 : ! Name of specs%method
734 0 : specs%method = 'Ionic positions relaxation using DIIS'
735 : ! Number of history
736 0 : specs%nhist = 3
737 : ! This is the initialization for ionmov==21
738 : ! -------------------------------------------
739 : case (21)
740 0 : specs%isARused=.TRUE.
741 : ! Values use in XML Output
742 0 : specs%type4xml='steepdesc'
743 0 : specs%crit4xml='tolmxf'
744 : ! Name of specs%method
745 0 : specs%method = 'Steepest descend algorithm'
746 : ! Number of history
747 0 : specs%nhist = 3
748 : ! This is the initialization for ionmov==22
749 : ! -------------------------------------------
750 : case (22)
751 : ! Values use in XML Output
752 2 : specs%type4xml='lbfgs'
753 2 : specs%crit4xml='tolmxf'
754 : ! Name of specs%method
755 2 : specs%method = 'Limited-memory Broyden-Fletcher-Goldfarb-Shanno method'
756 : ! Number of history
757 2 : specs%nhist = 3
758 : ! This is the initialization for ionmov==23
759 : ! -------------------------------------------
760 : case (23)
761 0 : specs%ncycle=2
762 : ! TEMPORARLY optcell is not allow
763 0 : specs%isARused=.FALSE.
764 0 : specs%isVused=.TRUE. ! Velocities are used
765 : ! Values use in XML Output
766 0 : specs%isFconv=.FALSE. ! Convergence is not used for MD
767 0 : specs%type4xml='isokin'
768 0 : specs%crit4xml='tolmxf'
769 : ! Name of specs%method
770 0 : specs%method = 'Using LOTF Molecular dynamics'
771 : ! Number of history
772 0 : specs%nhist = 3
773 : ! This is the initialization for ionmov==24
774 : ! -------------------------------------------
775 : case (24)
776 : specs%ncycle=1
777 : ! TEMPORARLY optcell is not allow
778 2 : specs%isARused=.FALSE.
779 2 : specs%isVused=.TRUE. ! Velocities are used
780 : ! Values use in XML Output
781 2 : specs%isFconv=.FALSE. ! Convergence is not used for MD
782 2 : specs%type4xml='velver'
783 2 : specs%crit4xml='none'
784 : ! Name of specs%method
785 2 : specs%method = 'Symplectic velocity verlet Molecular dynamics'
786 : ! Number of history
787 2 : specs%nhist = 3
788 : ! This is the initialization for ionmov==25
789 : ! -------------------------------------------
790 : case (25) ! Hybrid Monte Carlo algorithm (fixed lattice vectors)
791 2 : specs%ncycle = 12 ! Number of internal cycles (10+2)
792 2 : specs%isFconv=.FALSE. ! Convergence is not used for Monte Carlo
793 2 : specs%isVused=.TRUE. ! Velocities are used for update of atomic positions
794 : ! optcell is not allowed
795 2 : specs%isARused=.FALSE.
796 : ! Values use in XML Output
797 2 : specs%type4xml='hmc'
798 2 : specs%crit4xml='none'
799 : ! Name of specs%method
800 2 : specs%method = 'Hybrid Monte Carlo'
801 2 : specs%nhist = 3
802 : ! This is the initialization for ionmov==27
803 : ! -------------------------------------------
804 : case (27) ! Generation of the training set for effective potential
805 : specs%ncycle = 1 ! Number of internal cycles
806 0 : specs%isFconv=.FALSE. ! Convergence is not used
807 : specs%isVused=.FALSE. ! Velocities are not used for update of atomic positions
808 : ! Values use in XML Output
809 0 : specs%type4xml='TS'
810 0 : specs%crit4xml='none'
811 : ! Name of specs%method
812 0 : specs%method = 'training set generator'
813 : ! Number of history
814 0 : specs%nhist = -1
815 :
816 : case (28)
817 : ! Values used in XML Output
818 0 : specs%type4xml='simple'
819 0 : specs%crit4xml='tolmxf'
820 : ! Name of specs%method
821 0 : specs%method = "i-pi protocol"
822 : ! Number of history
823 : !specs%nhist = 3
824 : ! This is the initialization for ionmov==6
825 :
826 : case default
827 : ! MG TODO: Why is this check deactivated. We should have an empty case for all the ionmov that can use the default values
828 4368 : write(msg,"(a,i0)")"Wrong value for ionmov: ",ab_mover%ionmov
829 : !ABI_ERROR(msg)
830 : end select
831 :
832 4368 : end subroutine abimover_ini
833 : !!***
834 :
835 : !----------------------------------------------------------------------
836 :
837 : !!****f* m_abimover/abimover_destroy
838 : !! NAME
839 : !! abimover_destroy
840 : !!
841 : !! FUNCTION
842 : !! Destroy the abimover structure
843 : !!
844 : !! SIDE EFFECTS
845 : !! ab_mover <type(abimover)> = The abimover structure to be destroyed
846 : !!
847 : !! SOURCE
848 :
849 4368 : subroutine abimover_destroy(ab_mover)
850 :
851 : !Arguments ------------------------------------
852 : type(abimover),intent(inout) :: ab_mover
853 :
854 : ! ***************************************************************
855 :
856 4368 : nullify(ab_mover%goprecprm)
857 4368 : nullify(ab_mover%iatfix)
858 4368 : nullify(ab_mover%mdtemp)
859 4368 : nullify(ab_mover%ph_ngqpt)
860 4368 : nullify(ab_mover%ph_qshift)
861 :
862 4368 : nullify(ab_mover%prtatlist)
863 4368 : nullify(ab_mover%qmass)
864 4368 : nullify(ab_mover%strtarget)
865 4368 : nullify(ab_mover%symafm)
866 4368 : nullify(ab_mover%symrel)
867 4368 : nullify(ab_mover%tnons)
868 4368 : nullify(ab_mover%typat)
869 4368 : nullify(ab_mover%znucl)
870 :
871 4368 : nullify(ab_mover%amu_curr)
872 4368 : ABI_FREE(ab_mover%amass)
873 :
874 4368 : nullify(ab_mover%fnameabi_hes)
875 4368 : nullify(ab_mover%filnam_ds)
876 :
877 4368 : end subroutine abimover_destroy
878 : !!***
879 :
880 : !----------------------------------------------------------------------
881 :
882 : !!****f* defs_mover/abimover_print
883 : !! NAME
884 : !! abimover_print
885 : !!
886 : !! FUNCTION
887 : !! Print all the variables in a ab_mover
888 : !!
889 : !! INPUTS
890 : !!
891 : !! OUTPUT
892 : !!
893 : !! SIDE EFFECTS
894 : !! ab_mover <type(abimover)> = The ab_mover to nullify
895 : !!
896 : !! NOTES
897 : !! At present 29 variables are present in abimover
898 : !! if a new variable is added in abimover it should
899 : !! be added also for print here
900 : !!
901 : !! SOURCE
902 :
903 : subroutine abimover_print(ab_mover,iout)
904 :
905 : !Arguments ------------------------------------
906 : integer,intent(in) :: iout
907 : type(abimover),intent(inout) :: ab_mover
908 :
909 : !Local variables-------------------------------
910 : !arrays
911 : character(len=1200) :: message
912 : character(len=110) :: fmt
913 :
914 : ! ***************************************************************
915 :
916 : fmt='(a,e12.5,a,a,I5,a,a,I5,a,a,I5,a,a,I5,a,a,I5,a,a,I5,a,a,e12.5,a,a,e12.5,a,a,e12.5,a,a,e12.5,a,a,e12.5,a)'
917 :
918 : write(message,fmt)&
919 : & 'Delta Time for IONs',ab_mover%dtion,ch10, &
920 : & 'include a JELLium SLAB in the cell',ab_mover%jellslab,ch10, &
921 : & 'Number of ATOMs',ab_mover%natom,ch10, &
922 : & 'Number of CONstraint EQuations',ab_mover%nconeq,ch10, &
923 : & 'Number of SYMmetry operations',ab_mover%nsym,ch10, &
924 : & 'OPTimize the CELL shape and dimensions',ab_mover%optcell,ch10, &
925 : & 'RESTART Xcart and Fred',ab_mover%restartxf,ch10, &
926 : & 'Molecular Dynamics Initial Temperature',ab_mover%mdtemp(1),ch10, &
927 : & 'Molecular Dynamics Final Temperature',ab_mover%mdtemp(2),ch10, &
928 : & 'NOSE thermostat INERTia factor',ab_mover%noseinert,ch10, &
929 : & 'STRess PRECONditioner',ab_mover%strprecon,ch10, &
930 : & 'VIScosity',ab_mover%vis,ch10
931 :
932 : ! ! arrays
933 : ! ! Indices of AToms that are FIXed
934 : ! integer, pointer :: iatfix(:,:)
935 : ! ! SYMmetries, Anti-FerroMagnetic characteristics
936 : ! integer, pointer :: symafm(:)
937 : ! ! SYMmetry in REaL space
938 : ! integer, pointer :: symrel(:,:,:)
939 : ! Translation NON-Symmorphic vectors
940 : ! real(dp), pointer :: tnons(:,:)
941 : ! ! Mass of each atom (NOT IN DTSET)
942 : ! real(dp), pointer :: amass(:)
943 : ! ! STRess TARGET
944 : ! real(dp), pointer :: strtarget(:)
945 : ! Filename for Hessian matrix
946 : ! character(len=fnlen), pointer :: fnameabi_hes
947 :
948 : write(iout,*) 'CONTENT of ab_mover (scalar only)'
949 : write(iout,'(a)') message
950 :
951 : end subroutine abimover_print
952 : !!***
953 :
954 : !----------------------------------------------------------------------
955 :
956 : !!****f* m_abimover/mttk_ini
957 : !!
958 : !! NAME
959 : !! mttk_ini
960 : !!
961 : !! FUNCTION
962 : !! destructor function for mttk object
963 : !!
964 : !! INPUT
965 : !! mttk
966 : !!
967 : !! OUTPUT
968 : !!
969 : !! SOURCE
970 :
971 36 : subroutine mttk_ini(mttk_vars,nnos)
972 :
973 : integer,intent(in) :: nnos
974 : type(mttk_type), intent(out) :: mttk_vars
975 :
976 108 : ABI_MALLOC(mttk_vars%glogs,(nnos))
977 72 : ABI_MALLOC(mttk_vars%vlogs,(nnos))
978 72 : ABI_MALLOC(mttk_vars%xlogs,(nnos))
979 :
980 36 : end subroutine mttk_ini
981 : !!***
982 :
983 : !----------------------------------------------------------------------
984 :
985 : !!****f* m_abimover/mttk_fin
986 : !!
987 : !! NAME
988 : !! mttk_fin
989 : !!
990 : !! FUNCTION
991 : !! destructor function for mttk object
992 : !!
993 : !! INPUT
994 : !! mttk
995 : !!
996 : !! OUTPUT
997 : !!
998 : !! SOURCE
999 :
1000 4006 : subroutine mttk_fin(mttk_vars)
1001 :
1002 : type(mttk_type), intent(inout) :: mttk_vars
1003 :
1004 4006 : ABI_SFREE(mttk_vars%glogs)
1005 4006 : ABI_SFREE(mttk_vars%vlogs)
1006 4006 : ABI_SFREE(mttk_vars%xlogs)
1007 :
1008 4006 : end subroutine mttk_fin
1009 : !!***
1010 :
1011 : !----------------------------------------------------------------------
1012 :
1013 : !!****f* m_abimover/abiforstr_ini
1014 : !!
1015 : !! NAME
1016 : !! abiforstr_ini
1017 : !!
1018 : !! FUNCTION
1019 : !! destructor function for abiforstr object
1020 : !!
1021 : !! INPUT
1022 : !! forstr
1023 : !!
1024 : !! OUTPUT
1025 : !!
1026 : !! SOURCE
1027 :
1028 11919 : subroutine abiforstr_ini(forstr,natom)
1029 :
1030 : integer,intent(in) :: natom
1031 : type(abiforstr), intent(out) :: forstr
1032 :
1033 35757 : ABI_MALLOC(forstr%fcart,(3,natom))
1034 23838 : ABI_MALLOC(forstr%gred,(3,natom))
1035 :
1036 11919 : end subroutine abiforstr_ini
1037 : !!***
1038 :
1039 : !----------------------------------------------------------------------
1040 :
1041 : !!****f* m_abimover/abiforstr_fin
1042 : !!
1043 : !! NAME
1044 : !! abiforstr_fin
1045 : !!
1046 : !! FUNCTION
1047 : !! destructor function for abiforstr object
1048 : !!
1049 : !! INPUT
1050 : !! forstr
1051 : !!
1052 : !! OUTPUT
1053 : !!
1054 : !! SOURCE
1055 :
1056 11919 : subroutine abiforstr_fin(forstr)
1057 :
1058 : type(abiforstr), intent(inout) :: forstr
1059 :
1060 11919 : ABI_SFREE(forstr%fcart)
1061 11919 : ABI_SFREE(forstr%gred)
1062 :
1063 11919 : end subroutine abiforstr_fin
1064 : !!***
1065 :
1066 : !----------------------------------------------------------------------
1067 :
1068 : !!****f* m_abimover/make_prim_internals
1069 : !! NAME
1070 : !! make_prim_internals
1071 : !!
1072 : !! FUNCTION
1073 : !! Determine the bonds, angles and dihedrals for a starting
1074 : !! geometry, based on covalent radii for the atoms.
1075 : !!
1076 : !! INPUTS
1077 : !! natom = Number of atoms (dtset%natom)
1078 : !! nrshift= dimension of rshift
1079 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1080 : !! rshift(3,nrshift)=shift in xred that must be done to find all neighbors of
1081 : !! a given atom within a given number of neighboring shells
1082 : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
1083 : !!
1084 : !! OUTPUT
1085 : !!
1086 : !! SIDE EFFECTS
1087 : !! deloc <type(delocint)>=Important variables for
1088 : !! | pred_delocint
1089 : !! ! icenter = Index of the center of the number of shifts
1090 : !! | nang = Number of angles
1091 : !! | nbond = Number of bonds
1092 : !! | ncart = Number of cartesian directions
1093 : !! | (used for constraints)
1094 : !! | ndihed = Number of dihedrals
1095 : !! | nrshift = Dimension of rshift
1096 : !! | ninternal= Number of internal coordinates
1097 : !! | ninternal=nbond+nang+ndihed+ncart
1098 : !! |
1099 : !! | angs(2,3,nang) = Indexes to characterize angles
1100 : !! | bonds(2,2,nbond)= For a bond between iatom and jatom
1101 : !! | bonds(1,1,nbond) = iatom
1102 : !! | bonds(2,1,nbond) = icenter
1103 : !! | bonds(1,2,nbond) = jatom
1104 : !! | bonds(2,2,nbond) = irshift
1105 : !! | carts(2,ncart) = Index of total primitive internal,
1106 : !! | and atom (carts(2,:))
1107 : !! | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals
1108 : !! |
1109 : !! | rshift(3,nrshift)= Shift in xred that must be done to find
1110 : !! | all neighbors of a given atom within a
1111 : !! | given number of neighboring shells
1112 : !!
1113 : !! NOTES
1114 : !!
1115 : !! Adds cartesian coordinates if the number of internals with a
1116 : !! given atom is < 4 the chosen coordinate could be optimized
1117 : !! to be less dependent of the internals already incorporated.
1118 : !!
1119 : !! SOURCE
1120 :
1121 1 : subroutine make_prim_internals(deloc,natom,ntypat,rprimd,typat,xcart,znucl)
1122 :
1123 : !Arguments ------------------------------------
1124 : !scalars
1125 : type(delocint),intent(inout) :: deloc
1126 : integer,intent(in) :: natom,ntypat
1127 : !arrays
1128 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1129 : integer,intent(in) :: typat(natom)
1130 : real(dp),intent(in) :: znucl(:) ! znucl(ntypat) or znucl(npsp) ?
1131 :
1132 : !Local variables ------------------------------
1133 : ! function
1134 : !scalars
1135 : integer :: iang,iatom,ibond,icart,idihed,ii
1136 : real(dp) :: spp
1137 : !arrays
1138 2 : integer :: particip_atom(natom)
1139 1 : integer,allocatable :: badangles(:)
1140 : real(dp) :: rpt1(3),rpt3(3) !rpt2(3)
1141 :
1142 : !************************************************************************
1143 :
1144 3 : particip_atom(:) = 0
1145 :
1146 1 : call make_bonds(deloc,natom,ntypat,rprimd,typat,xcart,znucl)
1147 :
1148 5 : do ibond=1,deloc%nbond
1149 4 : write(std_out,'(a,i4,2(2i5,2x))') 'bond ', ibond, deloc%bonds(:,:,ibond)
1150 4 : particip_atom(deloc%bonds(1,1,ibond)) = particip_atom(deloc%bonds(1,1,ibond))+1
1151 5 : particip_atom(deloc%bonds(1,2,ibond)) = particip_atom(deloc%bonds(1,2,ibond))+1
1152 : end do
1153 :
1154 1 : call make_angles(deloc,natom)
1155 :
1156 3 : ABI_MALLOC(badangles,(deloc%nang))
1157 13 : badangles(:) = 0
1158 13 : do iang=1,deloc%nang
1159 12 : write(std_out,'(a,i4,3(2i5,2x))') 'angle ', iang, deloc%angs(:,:,iang)
1160 12 : particip_atom(deloc%angs(1,1,iang)) = particip_atom(deloc%angs(1,1,iang))+1
1161 12 : particip_atom(deloc%angs(1,2,iang)) = particip_atom(deloc%angs(1,2,iang))+1
1162 12 : particip_atom(deloc%angs(1,3,iang)) = particip_atom(deloc%angs(1,3,iang))+1
1163 :
1164 : ! DEBUG
1165 : ! rpt1(:) = xcart(:,deloc%angs(1,1,iang)) &
1166 : ! & + deloc%rshift(1,deloc%angs(2,1,iang))*rprimd(:,1) &
1167 : ! & + deloc%rshift(2,deloc%angs(2,1,iang))*rprimd(:,2) &
1168 : ! & + deloc%rshift(3,deloc%angs(2,1,iang))*rprimd(:,3)
1169 : ! rpt2(:) = xcart(:,deloc%angs(1,2,iang)) &
1170 : ! & + deloc%rshift(1,deloc%angs(2,2,iang))*rprimd(:,1) &
1171 : ! & + deloc%rshift(2,deloc%angs(2,2,iang))*rprimd(:,2) &
1172 : ! & + deloc%rshift(3,deloc%angs(2,2,iang))*rprimd(:,3)
1173 : ! rpt3(:) = xcart(:,deloc%angs(1,3,iang)) &
1174 : ! & + deloc%rshift(1,deloc%angs(2,3,iang))*rprimd(:,1) &
1175 : ! & + deloc%rshift(2,deloc%angs(2,3,iang))*rprimd(:,2) &
1176 : ! & + deloc%rshift(3,deloc%angs(2,3,iang))*rprimd(:,3)
1177 : ! write(std_out,*) rpt1,rpt2,rpt3,bond_length(rpt1,rpt2),bond_length(rpt2,rpt3)
1178 : ! ENDDEBUG
1179 :
1180 : ! check if angles are 180 degrees: discard the dihedrals in that case.
1181 : rpt1(:) = xcart(:,deloc%angs(1,1,iang)) &
1182 : & + deloc%rshift(1,deloc%angs(2,1,iang))*rprimd(:,1) &
1183 : & + deloc%rshift(2,deloc%angs(2,1,iang))*rprimd(:,2) &
1184 : & + deloc%rshift(3,deloc%angs(2,1,iang))*rprimd(:,3) &
1185 : & - xcart(:,deloc%angs(1,2,iang)) &
1186 : & - deloc%rshift(1,deloc%angs(2,2,iang))*rprimd(:,1) &
1187 : & - deloc%rshift(2,deloc%angs(2,2,iang))*rprimd(:,2) &
1188 48 : & - deloc%rshift(3,deloc%angs(2,2,iang))*rprimd(:,3)
1189 :
1190 : rpt3(:) = xcart(:,deloc%angs(1,3,iang)) &
1191 : & + deloc%rshift(1,deloc%angs(2,3,iang))*rprimd(:,1) &
1192 : & + deloc%rshift(2,deloc%angs(2,3,iang))*rprimd(:,2) &
1193 : & + deloc%rshift(3,deloc%angs(2,3,iang))*rprimd(:,3) &
1194 : & - xcart(:,deloc%angs(1,2,iang)) &
1195 : & - deloc%rshift(1,deloc%angs(2,2,iang))*rprimd(:,1) &
1196 : & - deloc%rshift(2,deloc%angs(2,2,iang))*rprimd(:,2) &
1197 48 : & - deloc%rshift(3,deloc%angs(2,2,iang))*rprimd(:,3)
1198 : spp = (rpt1(1)*rpt3(1)+rpt1(2)*rpt3(2)+rpt1(3)*rpt3(3))&
1199 : & / sqrt(rpt1(1)*rpt1(1)+rpt1(2)*rpt1(2)+rpt1(3)*rpt1(3)) &
1200 12 : & / sqrt(rpt3(1)*rpt3(1)+rpt3(2)*rpt3(2)+rpt3(3)*rpt3(3))
1201 13 : if (abs(abs(spp) - one) < tol6) then
1202 0 : write(std_out,*) 'make_prim_internals : an angle is too close to 180 degrees:'
1203 0 : write(std_out,*) ' will discard dihedrals using it '
1204 0 : badangles(iang) = 1
1205 : end if
1206 : end do
1207 :
1208 1 : call make_dihedrals(badangles,deloc)
1209 1 : ABI_FREE(badangles)
1210 :
1211 37 : do idihed=1,deloc%ndihed
1212 36 : write(std_out,'(a,i4,4(2i5,2x))') 'dihedral ', idihed, deloc%dihedrals(:,:,idihed)
1213 36 : particip_atom(deloc%dihedrals(1,1,idihed)) = particip_atom(deloc%dihedrals(1,1,idihed))+1
1214 36 : particip_atom(deloc%dihedrals(1,2,idihed)) = particip_atom(deloc%dihedrals(1,2,idihed))+1
1215 36 : particip_atom(deloc%dihedrals(1,3,idihed)) = particip_atom(deloc%dihedrals(1,3,idihed))+1
1216 37 : particip_atom(deloc%dihedrals(1,4,idihed)) = particip_atom(deloc%dihedrals(1,4,idihed))+1
1217 :
1218 : ! do ii=1,4
1219 : ! write(std_out,'((3E16.6,2x))') xcart(:,deloc%dihedrals(1,ii,idihed)) + &
1220 : ! & deloc%rshift(1,deloc%dihedrals(2,ii,idihed))*rprimd(:,1) + &
1221 : ! & deloc%rshift(2,deloc%dihedrals(2,ii,idihed))*rprimd(:,2) + &
1222 : ! & deloc%rshift(2,deloc%dihedrals(2,ii,idihed))*rprimd(:,3)
1223 : ! end do
1224 : end do
1225 :
1226 1 : write(std_out,*) 'make_deloc_internals: nbond,nang,ndihed = ', deloc%nbond,deloc%nang,deloc%ndihed
1227 :
1228 : !Check all atoms participate in at least 4 primitives. Otherwise, we should
1229 : !probably add cartesian coordinates to the internal ones.
1230 1 : deloc%ncart = 0
1231 3 : do iatom=1,natom
1232 3 : if (particip_atom(iatom) < 4) then
1233 0 : write(std_out,*) ' make_prim_internals : Warning : atom ', iatom, &
1234 0 : & ' does not belong to enough primitives to determine its'
1235 0 : write(std_out,*) ' position uniquely ! instead : ', particip_atom(iatom)
1236 0 : write(std_out,*) ' Will add cartesian coordinates to set of internals.'
1237 : ! write(std_out,*) ' Not done yet.'
1238 : ! stop
1239 0 : deloc%ncart = deloc%ncart + 4-particip_atom(iatom)
1240 : end if
1241 : end do
1242 1 : ABI_SFREE(deloc%carts)
1243 3 : ABI_MALLOC(deloc%carts ,(2,deloc%ncart))
1244 1 : icart = 0
1245 3 : do iatom=1,natom
1246 3 : if (particip_atom(iatom) < 4) then
1247 : ! kind of arbitrary : include first few directions for the atom: x, then y then z
1248 0 : do ii=1,4-particip_atom(iatom)
1249 0 : icart = icart+1
1250 0 : deloc%carts(1,icart) = ii
1251 0 : deloc%carts(2,icart) = iatom
1252 : end do
1253 : end if
1254 : end do
1255 :
1256 : !ninternal=nbond+nang+ndihed
1257 1 : deloc%ninternal=deloc%nbond+deloc%nang+deloc%ndihed+deloc%ncart
1258 :
1259 1 : end subroutine make_prim_internals
1260 : !!***
1261 :
1262 : !----------------------------------------------------------------------
1263 :
1264 : !!****f* m_abimover/make_angles
1265 : !! NAME
1266 : !! make_angles
1267 : !!
1268 : !! FUNCTION
1269 : !!
1270 : !! INPUTS
1271 : !!
1272 : !! OUTPUT
1273 : !!
1274 : !! SOURCE
1275 :
1276 1 : subroutine make_angles(deloc,natom)
1277 :
1278 : !Arguments ------------------------------------
1279 : !scalars
1280 : integer,intent(in) :: natom
1281 : type(delocint),intent(inout) :: deloc
1282 : !arrays
1283 :
1284 : !Local variables-------------------------------
1285 : !scalars
1286 : integer :: ia1,ia2,ibond,is1,is2,ishift,ja1,ja2
1287 : integer :: jbond,js1,js2
1288 : !arrays
1289 1 : integer,allocatable :: angs_tmp(:,:,:)
1290 :
1291 : ! *************************************************************************
1292 :
1293 : !tentative first allocation: < 6 angles per bond.
1294 3 : ABI_MALLOC(angs_tmp,(2,3,72*natom))
1295 :
1296 1 : deloc%nang = 0
1297 :
1298 5 : do ibond=1, deloc%nbond
1299 4 : ia1 = deloc%bonds(1,1,ibond)
1300 4 : is1 = deloc%bonds(2,1,ibond)
1301 4 : ia2 = deloc%bonds(1,2,ibond)
1302 4 : is2 = deloc%bonds(2,2,ibond)
1303 11 : do jbond=ibond+1,deloc%nbond
1304 6 : ja1 = deloc%bonds(1,1,jbond)
1305 6 : ja2 = deloc%bonds(1,2,jbond)
1306 2068 : do ishift=-(deloc%icenter-1),+(deloc%icenter-1)
1307 2058 : js1 = deloc%bonds(2,1,jbond)+ishift
1308 2058 : js2 = deloc%bonds(2,2,jbond)+ishift
1309 :
1310 2058 : if (ia1==ja1 .and. is1==js1) then
1311 6 : deloc%nang = deloc%nang+1
1312 18 : angs_tmp(:,1,deloc%nang) = (/ia2,is2/)
1313 18 : angs_tmp(:,2,deloc%nang) = (/ia1,is1/)
1314 18 : angs_tmp(:,3,deloc%nang) = (/ja2,js2/)
1315 :
1316 2052 : else if (ia1==ja2 .and. is1==js2) then
1317 0 : deloc%nang = deloc%nang+1
1318 0 : angs_tmp(:,1,deloc%nang) = (/ia2,is2/)
1319 0 : angs_tmp(:,2,deloc%nang) = (/ia1,is1/)
1320 0 : angs_tmp(:,3,deloc%nang) = (/ja1,js1/)
1321 :
1322 2052 : else if (ia2==ja2 .and. is2==js2) then
1323 6 : deloc%nang = deloc%nang+1
1324 18 : angs_tmp(:,1,deloc%nang) = (/ia1,is1/)
1325 18 : angs_tmp(:,2,deloc%nang) = (/ia2,is2/)
1326 18 : angs_tmp(:,3,deloc%nang) = (/ja1,js1/)
1327 :
1328 2046 : else if (ia2==ja1 .and. is2==js1) then
1329 0 : deloc%nang = deloc%nang+1
1330 0 : angs_tmp(:,1,deloc%nang) = (/ia1,is1/)
1331 0 : angs_tmp(:,2,deloc%nang) = (/ia2,is2/)
1332 0 : angs_tmp(:,3,deloc%nang) = (/ja2,js2/)
1333 :
1334 : end if
1335 2064 : if (deloc%nang > 72*natom) then
1336 0 : ABI_ERROR('too many angles found > 72*natom')
1337 : end if
1338 : end do
1339 : end do ! jbond do
1340 : end do ! ibond
1341 :
1342 1 : ABI_SFREE(deloc%angs)
1343 3 : ABI_MALLOC(deloc%angs,(2,3,deloc%nang))
1344 121 : deloc%angs(:,:,1:deloc%nang) = angs_tmp(:,:,1:deloc%nang)
1345 1 : ABI_FREE(angs_tmp)
1346 :
1347 1 : end subroutine make_angles
1348 : !!***
1349 :
1350 : !----------------------------------------------------------------------
1351 :
1352 : !!****f* m_abimover/make_dihedrals
1353 : !! NAME
1354 : !! make_dihedrals
1355 : !!
1356 : !! FUNCTION
1357 : !!
1358 : !! INPUTS
1359 : !!
1360 : !! OUTPUT
1361 : !!
1362 : !! SOURCE
1363 :
1364 1 : subroutine make_dihedrals(badangles,deloc)
1365 :
1366 : !Arguments ------------------------------------
1367 : !scalars
1368 : type(delocint),intent(inout) :: deloc
1369 : !arrays
1370 : integer,intent(in) :: badangles(deloc%nang)
1371 :
1372 : !Local variables-------------------------------
1373 : !scalars
1374 : integer :: chkdihed,ia1,ia2,ia3,iang,idihed,is1,is2
1375 : integer :: is3,ishift,ja1,ja2,ja3,jang,js1,js2,js3,maxshift
1376 : integer :: minshift
1377 : !arrays
1378 1 : integer,allocatable :: diheds_tmp(:,:,:)
1379 :
1380 : ! *************************************************************************
1381 :
1382 : !tentative first allocation: < 6 dihedrals per angle.
1383 3 : ABI_MALLOC(diheds_tmp,(2,4,6*deloc%nang))
1384 :
1385 1 : deloc%ndihed = 0
1386 937 : diheds_tmp(:,:,:) = 0
1387 :
1388 13 : do iang=1,deloc%nang
1389 12 : if (badangles(iang) == 1) cycle
1390 12 : ia1 = deloc%angs(1,1,iang)
1391 12 : is1 = deloc%angs(2,1,iang)
1392 12 : ia2 = deloc%angs(1,2,iang)
1393 12 : is2 = deloc%angs(2,2,iang)
1394 12 : ia3 = deloc%angs(1,3,iang)
1395 12 : is3 = deloc%angs(2,3,iang)
1396 :
1397 79 : do jang=iang+1,deloc%nang
1398 66 : if (badangles(jang) == 1) cycle
1399 66 : ja1 = deloc%angs(1,1,jang)
1400 66 : ja2 = deloc%angs(1,2,jang)
1401 66 : ja3 = deloc%angs(1,3,jang)
1402 22716 : do ishift=-(deloc%icenter-1),(deloc%icenter-1)
1403 22638 : js1 = deloc%angs(2,1,jang)+ishift
1404 22638 : js2 = deloc%angs(2,2,jang)+ishift
1405 22638 : js3 = deloc%angs(2,3,jang)+ishift
1406 :
1407 22638 : chkdihed=0
1408 22638 : if (ia2==ja1 .and. is2==js1) then
1409 36 : if (ia1==ja2 .and. is1==js2) then
1410 14 : deloc%ndihed = deloc%ndihed+1
1411 42 : diheds_tmp(:,1,deloc%ndihed) = (/ia3,is3/)
1412 42 : diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/)
1413 42 : diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/)
1414 42 : diheds_tmp(:,4,deloc%ndihed) = (/ja3,js3/)
1415 14 : chkdihed=1
1416 22 : else if (ia3==ja2 .and. is3==js2) then
1417 8 : deloc%ndihed = deloc%ndihed+1
1418 24 : diheds_tmp(:,1,deloc%ndihed) = (/ia1,is1/)
1419 24 : diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/)
1420 24 : diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/)
1421 24 : diheds_tmp(:,4,deloc%ndihed) = (/ja3,js3/)
1422 8 : chkdihed=1
1423 : end if
1424 22602 : else if (ia2==ja3 .and. is2==js3) then
1425 36 : if (ia1==ja2 .and. is1==js2) then
1426 0 : deloc%ndihed = deloc%ndihed+1
1427 0 : diheds_tmp(:,1,deloc%ndihed) = (/ia3,is3/)
1428 0 : diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/)
1429 0 : diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/)
1430 0 : diheds_tmp(:,4,deloc%ndihed) = (/ja1,js1/)
1431 0 : chkdihed=1
1432 36 : else if (ia3==ja2 .and. is3==js2) then
1433 14 : deloc%ndihed = deloc%ndihed+1
1434 42 : diheds_tmp(:,1,deloc%ndihed) = (/ia1,is1/)
1435 42 : diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/)
1436 42 : diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/)
1437 42 : diheds_tmp(:,4,deloc%ndihed) = (/ja1,js1/)
1438 14 : chkdihed=1
1439 : end if
1440 : end if
1441 22638 : if (deloc%ndihed > 6*deloc%nang) then
1442 0 : ABI_ERROR('make_dihedrals : too many dihedrals found > 6*nang')
1443 : end if
1444 22704 : if (chkdihed == 1) then
1445 36 : if ( diheds_tmp(1,4,deloc%ndihed) == diheds_tmp(1,1,deloc%ndihed) .and.&
1446 : & diheds_tmp(2,4,deloc%ndihed) == diheds_tmp(2,1,deloc%ndihed) ) then
1447 0 : write(std_out,*) 'make_dihedrals : Bad dihedral was found: atom1 == atom4. Discarding.'
1448 0 : diheds_tmp(:,:,deloc%ndihed) = 0
1449 0 : deloc%ndihed = deloc%ndihed-1
1450 : end if
1451 : end if
1452 : end do
1453 : end do
1454 : ! end jang do
1455 : end do
1456 : !end iang do
1457 :
1458 1 : ABI_SFREE(deloc%dihedrals)
1459 :
1460 3 : ABI_MALLOC(deloc%dihedrals,(2,4,deloc%ndihed))
1461 37 : do idihed=1,deloc%ndihed
1462 468 : deloc%dihedrals(:,:,idihed) = diheds_tmp(:,:,idihed)
1463 :
1464 : ! minshift = minval(diheds_tmp(2,:,idihed))
1465 : ! if (minshift <= 0) then
1466 : ! deloc%dihedrals(2,:,idihed) = deloc%dihedrals(2,:,idihed)+minshift+1
1467 : ! end if
1468 : ! maxshift = maxval(diheds_tmp(2,:,idihed))
1469 : ! if (maxshift > deloc%nrshift) then
1470 : ! deloc%dihedrals(2,:,idihed) = deloc%dihedrals(2,:,idihed)-maxshift
1471 : ! end if
1472 : !
1473 180 : minshift = minval(diheds_tmp(2,:,idihed))
1474 180 : maxshift = maxval(diheds_tmp(2,:,idihed))
1475 37 : if (minshift <= 0 .or. maxshift > deloc%nrshift) then
1476 0 : ABI_ERROR("dihedral extends beyond first neighboring unit cells!")
1477 : end if
1478 : end do
1479 1 : ABI_FREE(diheds_tmp)
1480 :
1481 1 : end subroutine make_dihedrals
1482 : !!***
1483 :
1484 : !----------------------------------------------------------------------
1485 :
1486 : !!****f* m_abimover/make_bonds
1487 : !! NAME
1488 : !! make_bonds
1489 : !!
1490 : !! FUNCTION
1491 : !!
1492 : !! INPUTS
1493 : !!
1494 : !! OUTPUT
1495 : !!
1496 : !! SOURCE
1497 :
1498 1 : subroutine make_bonds(deloc,natom,ntypat,rprimd,typat,xcart,znucl)
1499 :
1500 : !Arguments ------------------------------------
1501 : !scalars
1502 : integer,intent(in) :: natom,ntypat
1503 : type(delocint),intent(inout) :: deloc
1504 : !arrays
1505 : integer,intent(in) :: typat(natom)
1506 : real(dp),intent(in) :: znucl(:) ! znucl(ntypat) or
1507 : ! znucl(npsp) ?
1508 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1509 :
1510 : !Local variables ------------------------------
1511 : !scalars
1512 : integer :: iatom,irshift,itypat,jatom
1513 : real(dp) :: bl,bondfudge,rcov1,rcov2
1514 : type(atomdata_t) :: atom
1515 : !arrays
1516 1 : integer,allocatable :: bonds_tmp(:,:,:)
1517 2 : real(dp) :: rcov(ntypat),rpt(3)
1518 :
1519 : !************************************************************************
1520 :
1521 2 : do itypat=1,ntypat
1522 1 : call atomdata_from_znucl(atom,znucl(itypat))
1523 2 : rcov(itypat) = atom%rcov
1524 : end do
1525 :
1526 : !write(std_out,*) ' rcov =', rcov
1527 : !write(std_out,*) ' nrshift =', deloc%nrshift
1528 : !write(std_out,*) ' xcart =', xcart
1529 : !write(std_out,*) ' natom =',natom
1530 :
1531 : !tentative first allocation: < 12 bonds per atom.
1532 3 : ABI_MALLOC(bonds_tmp,(2,2,12*natom))
1533 :
1534 1 : bondfudge = 1.1_dp
1535 :
1536 1 : deloc%nbond = 0
1537 :
1538 3 : do iatom=1,natom
1539 2 : rcov1 = rcov(typat(iatom))
1540 4 : do jatom=iatom+1,natom
1541 1 : rcov2 = rcov(typat(jatom))
1542 346 : do irshift=1,deloc%nrshift
1543 : rpt(:) = xcart(:,jatom) &
1544 : & + deloc%rshift(1,irshift)*rprimd(:,1) &
1545 : & + deloc%rshift(2,irshift)*rprimd(:,2) &
1546 1372 : & + deloc%rshift(3,irshift)*rprimd(:,3)
1547 343 : bl = bond_length(xcart(:,iatom),rpt)
1548 :
1549 : !write(std_out,*) ' bl, bondfudge*(rcov1+rcov2) = ',bl, bondfudge*(rcov1+rcov2)
1550 :
1551 344 : if (bondfudge*(rcov1+rcov2) - bl > tol6) then
1552 4 : deloc%nbond = deloc%nbond+1
1553 4 : if (deloc%nbond > 12*natom) then
1554 0 : ABI_ERROR('make_bonds: error too many bonds !')
1555 : end if
1556 4 : bonds_tmp(1,1,deloc%nbond) = iatom
1557 4 : bonds_tmp(2,1,deloc%nbond) = deloc%icenter
1558 4 : bonds_tmp(1,2,deloc%nbond) = jatom
1559 4 : bonds_tmp(2,2,deloc%nbond) = irshift
1560 :
1561 : !write(std_out,*) ' ibond bonds = ', deloc%nbond, bonds_tmp(:,:,deloc%nbond),xcart(:,iatom),rpt
1562 : end if
1563 : end do ! jatom
1564 : end do
1565 : end do ! iatom
1566 :
1567 1 : ABI_SFREE(deloc%bonds)
1568 :
1569 3 : ABI_MALLOC(deloc%bonds,(2,2,deloc%nbond))
1570 29 : deloc%bonds(:,:,1:deloc%nbond) = bonds_tmp(:,:,1:deloc%nbond)
1571 :
1572 : ! do ibond=1,deloc%nbond
1573 : ! write(std_out,*) ' make_bonds : bonds_tmp ', ibond, bonds_tmp(:,:,ibond)
1574 : ! write(std_out,*) ' make_bonds : deloc%bonds ', ibond, deloc%bonds(:,:,ibond)
1575 : ! end do
1576 :
1577 1 : ABI_FREE(bonds_tmp)
1578 :
1579 1 : end subroutine make_bonds
1580 : !!***
1581 :
1582 : !----------------------------------------------------------------------
1583 :
1584 : !!****f* m_abimover/calc_prim_int
1585 : !! NAME
1586 : !! calc_prim_int
1587 : !!
1588 : !! FUNCTION
1589 : !! calculate values of primitive internal coordinates as a function of
1590 : !! cartesian ones.
1591 : !!
1592 : !! INPUTS
1593 : !! angs= number of angles
1594 : !! bonds(2,2,nbond)=for a bond between iatom and jatom
1595 : !! bonds(1,1,nbond) = iatom
1596 : !! bonds(2,1,nbond) = icenter
1597 : !! bonds(1,2,nbond) = jatom
1598 : !! bonds(2,2,nbond) = irshift
1599 : !! carts(2,ncart)= index of total primitive internal, and atom (carts(2,:))
1600 : !! dihedrals(2,4,ndihed)=indexes to characterize dihedrals
1601 : !! dtset <type(dataset_type)>=all input variables for this dataset
1602 : !! nang(2,3,nang)=indexes to characterize angles
1603 : !! nbond=number of bonds
1604 : !! ncart=number of cartesian coordinates used
1605 : !! ndihed= number of dihedrals
1606 : !! ninternal=nbond+nang+ndihed+ncart: number of internal coordinates
1607 : !! nrshift= dimension of rshift
1608 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1609 : !! rshift(3,nrshift)=shift in xred that must be done to find all neighbors of
1610 : !! a given atom within a given number of neighboring shells
1611 : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
1612 : !!
1613 : !! OUTPUT
1614 : !! prim_int(ninternal)=values of primitive internal coordinates
1615 : !!
1616 : !! SIDE EFFECTS
1617 : !!
1618 : !! NOTES
1619 : !!
1620 : !! SOURCE
1621 :
1622 65 : subroutine calc_prim_int(deloc,natom,rprimd,xcart,prim_int)
1623 :
1624 : !Arguments ------------------------------------
1625 : !scalars
1626 : integer,intent(in) :: natom
1627 : type(delocint),intent(in) :: deloc
1628 : !arrays
1629 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1630 : real(dp),intent(out) :: prim_int(deloc%ninternal)
1631 :
1632 : !Local variables ------------------------------
1633 : !scalars
1634 : integer :: i1,i2,i3,i4,iang,ibond,icart,idihed,iprim,s1,s2,s3,s4
1635 : !arrays
1636 : real(dp) :: r1(3),r2(3),r3(3),r4(3)
1637 :
1638 : !************************************************************************
1639 :
1640 : !DEBUG
1641 : !write(std_out,*) ' calc_prim_int : enter'
1642 : !write(std_out,*) shape(deloc%bonds)
1643 : !do ibond=1,deloc%nbond
1644 : !do i1=1,2
1645 : !write(std_out,'(2I5)') deloc%bonds(:,i1,ibond)
1646 : !end do
1647 : !end do
1648 : !ENDDEBUG
1649 :
1650 65 : iprim=1
1651 : !first: bond values
1652 325 : do ibond=1,deloc%nbond
1653 260 : i1 = deloc%bonds(1,1,ibond)
1654 260 : s1 = deloc%bonds(2,1,ibond)
1655 : r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
1656 : & +deloc%rshift(2,s1)*rprimd(:,2)&
1657 1040 : & +deloc%rshift(3,s1)*rprimd(:,3)
1658 260 : i2 = deloc%bonds(1,2,ibond)
1659 260 : s2 = deloc%bonds(2,2,ibond)
1660 : r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
1661 : & +deloc%rshift(2,s2)*rprimd(:,2)&
1662 1040 : & +deloc%rshift(3,s2)*rprimd(:,3)
1663 260 : prim_int(iprim) = bond_length(r1,r2)
1664 325 : iprim=iprim+1
1665 : end do
1666 :
1667 : !second: angle values (ang)
1668 845 : do iang=1,deloc%nang
1669 780 : i1 = deloc%angs(1,1,iang)
1670 780 : s1 = deloc%angs(2,1,iang)
1671 : r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
1672 : & +deloc%rshift(2,s1)*rprimd(:,2)&
1673 3120 : & +deloc%rshift(3,s1)*rprimd(:,3)
1674 780 : i2 = deloc%angs(1,2,iang)
1675 780 : s2 = deloc%angs(2,2,iang)
1676 : r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
1677 : & +deloc%rshift(2,s2)*rprimd(:,2)&
1678 3120 : & +deloc%rshift(3,s2)*rprimd(:,3)
1679 780 : i3 = deloc%angs(1,3,iang)
1680 780 : s3 = deloc%angs(2,3,iang)
1681 : r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
1682 : & +deloc%rshift(2,s3)*rprimd(:,2)&
1683 3120 : & +deloc%rshift(3,s3)*rprimd(:,3)
1684 780 : prim_int(iprim) = angle_ang(r1,r2,r3)
1685 845 : iprim=iprim+1
1686 : end do
1687 :
1688 : !third: dihedral values
1689 2405 : do idihed=1,deloc%ndihed
1690 2340 : i1 = deloc%dihedrals(1,1,idihed)
1691 2340 : s1 = deloc%dihedrals(2,1,idihed)
1692 : r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
1693 : & +deloc%rshift(2,s1)*rprimd(:,2)&
1694 9360 : & +deloc%rshift(3,s1)*rprimd(:,3)
1695 2340 : i2 = deloc%dihedrals(1,2,idihed)
1696 2340 : s2 = deloc%dihedrals(2,2,idihed)
1697 : r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
1698 : & +deloc%rshift(2,s2)*rprimd(:,2)&
1699 9360 : & +deloc%rshift(3,s2)*rprimd(:,3)
1700 2340 : i3 = deloc%dihedrals(1,3,idihed)
1701 2340 : s3 = deloc%dihedrals(2,3,idihed)
1702 : r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
1703 : & +deloc%rshift(2,s3)*rprimd(:,2)&
1704 9360 : & +deloc%rshift(3,s3)*rprimd(:,3)
1705 2340 : i4 = deloc%dihedrals(1,4,idihed)
1706 2340 : s4 = deloc%dihedrals(2,4,idihed)
1707 : r4(:) = xcart(:,i4)+deloc%rshift(1,s4)*rprimd(:,1)&
1708 : & +deloc%rshift(2,s4)*rprimd(:,2)&
1709 9360 : & +deloc%rshift(3,s4)*rprimd(:,3)
1710 2340 : prim_int(iprim) = angle_dihedral(r1,r2,r3,r4)
1711 2405 : iprim=iprim+1
1712 : end do
1713 :
1714 65 : do icart=1,deloc%ncart
1715 0 : prim_int(iprim) = xcart(deloc%carts(1,icart),deloc%carts(2,icart))
1716 65 : iprim=iprim+1
1717 : end do
1718 :
1719 : !DEBUG
1720 : !write(std_out,*) 'Primitive internal coordinate values:'
1721 : !do iprim=1,ninternal
1722 : !if (iprim <= deloc%nbond) then
1723 : !write(std_out,*) iprim, prim_int(iprim)
1724 : !else
1725 : !write(std_out,*) iprim, prim_int(iprim), prim_int(iprim)/pi*180.0_dp
1726 : !end if
1727 : !end do
1728 : !ENDDEBUG
1729 :
1730 65 : end subroutine calc_prim_int
1731 : !!***
1732 :
1733 : !----------------------------------------------------------------------
1734 :
1735 : !!****f* m_abimover/bond_length
1736 : !! NAME
1737 : !! bond_length
1738 : !!
1739 : !! FUNCTION
1740 : !!
1741 : !! INPUTS
1742 : !!
1743 : !! OUTPUT
1744 : !!
1745 : !! SOURCE
1746 :
1747 7601 : pure function bond_length(r1,r2)
1748 :
1749 : !Arguments ------------------------------------
1750 : !scalars
1751 : real(dp) :: bond_length
1752 : !arrays
1753 : real(dp),intent(in) :: r1(3),r2(3)
1754 :
1755 : !Local variables ------------------------------------
1756 : !arrays
1757 : real(dp) :: rpt(3)
1758 :
1759 : !******************************************************************
1760 30404 : rpt(:) = r1(:)-r2(:)
1761 30404 : bond_length = norm2(rpt)
1762 :
1763 7601 : end function bond_length
1764 : !!***
1765 :
1766 : !!****f* m_abimover/angle_ang
1767 : !! NAME
1768 : !! angle_ang
1769 : !!
1770 : !! FUNCTION
1771 : !!
1772 : !! INPUTS
1773 : !!
1774 : !! OUTPUT
1775 : !!
1776 : !! SOURCE
1777 :
1778 780 : pure function angle_ang(r1,r2,r3)
1779 :
1780 : !Arguments ------------------------------------
1781 : !scalars
1782 : real(dp) :: angle_ang
1783 : !arrays
1784 : real(dp),intent(in) :: r1(3),r2(3),r3(3)
1785 :
1786 : !Local variables ------------------------------
1787 : !scalars
1788 : real(dp) :: cos_ang,n1,n2
1789 : !arrays
1790 : real(dp) :: rpt12(3),rpt32(3)
1791 :
1792 : !******************************************************************
1793 780 : n1=bond_length(r1,r2)
1794 780 : n2=bond_length(r3,r2)
1795 :
1796 3120 : rpt12(:) = r1(:)-r2(:)
1797 3120 : rpt32(:) = r3(:)-r2(:)
1798 :
1799 3120 : cos_ang = dot_product(rpt12, rpt32) / (n1 * n2)
1800 :
1801 780 : if (cos_ang > one - epsilon(one)*two) then
1802 : cos_ang = one
1803 780 : else if(cos_ang < -one + epsilon(one)*two) then
1804 0 : cos_ang = -one
1805 : end if
1806 :
1807 780 : angle_ang=acos(cos_ang)
1808 :
1809 780 : end function angle_ang
1810 : !!***
1811 :
1812 : !!****f* m_abimover/angle_dihedral
1813 : !! NAME
1814 : !! angle_dihedral
1815 : !!
1816 : !! FUNCTION
1817 : !!
1818 : !! INPUTS
1819 : !!
1820 : !! OUTPUT
1821 : !!
1822 : !! SOURCE
1823 :
1824 2340 : function angle_dihedral(r1,r2,r3,r4)
1825 :
1826 : !Arguments ------------------------------------
1827 : !scalars
1828 : real(dp) :: angle_dihedral
1829 : !arrays
1830 : real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
1831 :
1832 : !Local variables------------------------------------
1833 : !scalars
1834 : real(dp) :: cos_dihedral,dih_sign,n1,n2,sin_dihedral
1835 : !arrays
1836 : real(dp) :: cp1232(3),cp3432(3),cpcp(3),rpt12(3),rpt32(3),rpt34(3)
1837 :
1838 : !******************************************************************
1839 :
1840 9360 : rpt12(:) = r1(:)-r2(:)
1841 9360 : rpt32(:) = r3(:)-r2(:)
1842 9360 : rpt34(:) = r3(:)-r4(:)
1843 :
1844 2340 : call acrossb(rpt12,rpt32,cp1232)
1845 2340 : call acrossb(rpt34,rpt32,cp3432)
1846 :
1847 : !DEBUG
1848 : !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
1849 : !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
1850 : !ENDDEBUG
1851 :
1852 9360 : n1 = norm2(cp1232)
1853 9360 : n2 = norm2(cp3432)
1854 :
1855 9360 : cos_dihedral = dot_product(cp1232, cp3432) / (n1 * n2)
1856 : !we use complementary of standard angle, so
1857 2340 : cos_dihedral = -cos_dihedral
1858 :
1859 2340 : call acrossb(cp1232,cp3432,cpcp)
1860 9360 : cpcp(:) = cpcp(:)/n1/n2
1861 16380 : sin_dihedral = -dot_product(cpcp, rpt32) / norm2(rpt32)
1862 2340 : dih_sign = one
1863 : !if (abs(sin_dihedral) > tol12) then
1864 : !dih_sign = sin_dihedral/abs(sin_dihedral)
1865 : !end if
1866 2340 : if (sin_dihedral < -tol12) then
1867 780 : dih_sign = -one
1868 : end if
1869 :
1870 : !DEBUG
1871 : !write(std_out,'(a,3E20.10)') 'angle_dihedral : cos sin dih_sign= ',&
1872 : !& cos_dihedral,sin_dihedral,dih_sign
1873 : !ENDDEBUG
1874 :
1875 2340 : if (cos_dihedral > one - epsilon(one)*two) then
1876 : cos_dihedral = one
1877 2340 : else if(cos_dihedral < -one + epsilon(one)*two) then
1878 774 : cos_dihedral = -one
1879 : end if
1880 :
1881 2340 : angle_dihedral = dih_sign*acos(cos_dihedral)
1882 :
1883 2340 : end function angle_dihedral
1884 : !!***
1885 :
1886 : !----------------------------------------------------------------------
1887 :
1888 : !!****f* m_abimover/make_bonds_new
1889 : !! NAME
1890 : !! make_bonds_new
1891 : !!
1892 : !! FUNCTION
1893 : !! Fill the contents of the bonds structure, that contains
1894 : !! all non redundant bonds that could be generated between
1895 : !! all the atoms in the unitary cell and their adjacent cells
1896 : !!
1897 : !! INPUTS
1898 : !! natom= Number of atoms
1899 : !! ntypat= Number of type of atoms
1900 : !! rprimd= Dimensional primitive vectors of the cell
1901 : !! xcart= Cartesian coordinates of the atoms
1902 : !! znucl= Z number of the atom
1903 : !!
1904 : !! OUTPUT
1905 : !! bonds= Structure that store all the information about
1906 : !! bonds created by this routine:
1907 : !! nbonds= Total number of bonds
1908 : !! nbondi= Number of bonds for atom i
1909 : !! indexi= Indeces of bonds for atom i
1910 : !! bond_length= Distances between atoms i and j (including shift)
1911 : !! bond_vect= Unitary vector for the bond from i to j
1912 : !! tolerance= The tolerance is multiplied to the
1913 : !! adition of covalent radius to decide if a bond is created
1914 : !!
1915 : !! SOURCE
1916 :
1917 6 : subroutine make_bonds_new(bonds,natom,ntypat,rprimd,typat,xcart,znucl)
1918 :
1919 : !Arguments ------------------------------------
1920 : !scalars
1921 : integer,intent(in) :: natom,ntypat
1922 : !arrays
1923 : integer,intent(in) :: typat(natom)
1924 : real(dp),intent(in) :: znucl(ntypat)
1925 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1926 : type(go_bonds),intent(inout) :: bonds
1927 :
1928 : !Local variables ------------------------------
1929 : !scalars
1930 : integer :: ii,jj,kk,ibond,irshift
1931 : real(dp) :: rcov1,rcov2
1932 : real(dp) :: bl
1933 6 : type(go_bonds) :: bonds_tmp
1934 : type(atomdata_t) :: atom
1935 :
1936 : !arrays
1937 12 : character(len=2) :: symbol(ntypat)
1938 12 : real(dp) :: amu(ntypat)
1939 : integer :: shift(3,13) ! Represent all shift vectors that are not equivalent by central symmetry
1940 : ! For example (1,1,-1) is equivalent to (-1,-1,1)
1941 : ! It means that bond between atom i in the original cell and atom j in the
1942 : ! cell with cordinates (1,1,-1) is equivalent to the bond between atom j in
1943 : ! the orignal cell and atom i in the cell with coordinates (-1,-1,1)
1944 : ! The trivial shift (0,0,0) is excluded here
1945 12 : real(dp) :: rcov(ntypat) ! Covalent radius
1946 : real(dp) :: rpt(3)
1947 :
1948 : !***************************************************************************
1949 : !Beginning of executable session
1950 : !***************************************************************************
1951 :
1952 : !write(std_out,*) 'make_bonds 01'
1953 : !##########################################################
1954 : !### 01. Compute covalent radius
1955 :
1956 12 : do ii=1,ntypat
1957 6 : call atomdata_from_znucl(atom,znucl(ii))
1958 6 : amu(ii) = atom%amu
1959 6 : rcov(ii) = atom%rcov
1960 12 : symbol(ii) = atom%symbol
1961 : end do
1962 :
1963 : !write(std_out,*) 'make_bonds 02'
1964 : !##########################################################
1965 : !### 02. Fill the 13 posible shift conecting adjacent cells
1966 :
1967 : shift(:,:)=reshape( (/ 1,0,0,&
1968 : & 0, 1, 0,&
1969 : & 0, 0, 1,&
1970 : & 1, 1, 0,&
1971 : & 1,-1, 0,&
1972 : & 0, 1, 1,&
1973 : & 0, 1,-1,&
1974 : & 1, 0, 1,&
1975 : & 1, 0,-1,&
1976 : & 1, 1, 1,&
1977 : & 1,-1, 1,&
1978 : & 1, 1,-1,&
1979 6 : & 1,-1,-1 /), (/ 3, 13 /))
1980 :
1981 : !write(std_out,*) 'make_bonds 03'
1982 : !##########################################################
1983 : !### 03. Initialize the values of bonds
1984 :
1985 : !The total number of bonds could not be predicted without
1986 : !compute all the distances, but the extreme case is linking
1987 : !all the atoms within all adjacent cells (natom*natom*13)
1988 : !plus the all the bonds inside the original cell (natom*(natom-1))
1989 :
1990 6 : bonds_tmp%nbonds=0
1991 6 : bonds_tmp%tolerance=bonds%tolerance
1992 6 : ibond=0
1993 :
1994 18 : ABI_MALLOC(bonds_tmp%bond_vect,(3,natom*natom*14-natom))
1995 18 : ABI_MALLOC(bonds_tmp%bond_length,(natom*natom*14-natom))
1996 :
1997 : !indexi contains the indices to the bonds
1998 24 : ABI_MALLOC(bonds_tmp%indexi,(natom,natom*natom*14-natom))
1999 :
2000 18 : ABI_MALLOC(bonds_tmp%nbondi,(natom))
2001 :
2002 978 : bonds_tmp%indexi(:,:)=0
2003 18 : bonds_tmp%nbondi(:)=0
2004 1302 : bonds_tmp%bond_vect(:,:)=0.0
2005 330 : bonds_tmp%bond_length(:)=0.0
2006 :
2007 : !write(std_out,*) 'make_bonds 04'
2008 : !##########################################################
2009 : !### 04. Compute the bonds inside the original cell
2010 : !### shift=(0,0,0)
2011 :
2012 18 : do ii=1,natom
2013 12 : rcov1 = rcov(typat(ii))
2014 :
2015 24 : do jj=ii+1,natom
2016 6 : rcov2 = rcov(typat(jj))
2017 :
2018 6 : bl=bond_length(xcart(:,ii),xcart(:,jj))
2019 :
2020 18 : if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then
2021 : ! We have a new bond, nbonds starts from
2022 : ! 0, so it could be used to index the
2023 : ! locations of bondij and distij
2024 :
2025 : ! Increase the number of bonds
2026 6 : bonds_tmp%nbonds= bonds_tmp%nbonds+1
2027 :
2028 : ! The number of bonds for atoms ii and jj
2029 : ! needs to raise by one
2030 6 : bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1
2031 6 : bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1
2032 :
2033 6 : bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds
2034 : ! The value for jj is negative to indicate that
2035 : ! the vector is from ii to jj
2036 6 : bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds
2037 :
2038 : ! The unitary vector is always from ii to jj
2039 24 : bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(xcart(:,jj)-xcart(:,ii))/bl
2040 6 : bonds_tmp%bond_length(bonds_tmp%nbonds)=bl
2041 :
2042 : end if
2043 :
2044 : end do !! jj
2045 : end do !! ii
2046 :
2047 : !write(std_out,*) 'make_bonds 05'
2048 : !##########################################################
2049 : !### 05. Compute the bonds outside the original cell
2050 : !### 13 shifts considered
2051 :
2052 : !Bonds between identical atoms but in diferent cells are
2053 : !allowed
2054 :
2055 18 : do ii=1,natom
2056 12 : rcov1 = rcov(typat(ii))
2057 42 : do jj=1,natom
2058 24 : rcov2 = rcov(typat(jj))
2059 :
2060 348 : do irshift=1,13
2061 :
2062 1248 : do kk=1,3
2063 : rpt(kk) = xcart(kk,jj)+&
2064 : & shift(1,irshift)*rprimd(kk,1)+ &
2065 : & shift(2,irshift)*rprimd(kk,2)+ &
2066 1248 : & shift(3,irshift)*rprimd(kk,3)
2067 : end do
2068 :
2069 :
2070 312 : bl =bond_length(xcart(:,ii),rpt)
2071 :
2072 336 : if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then
2073 :
2074 : ! We have a new bond, nbonds starts from
2075 : ! 0, so it could be used to index the
2076 : ! locations of bondij and distij
2077 :
2078 : ! Increase the number of bonds
2079 0 : bonds_tmp%nbonds= bonds_tmp%nbonds+1
2080 :
2081 : ! The number of bonds for atoms ii and jj
2082 : ! needs to raise by one
2083 0 : bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1
2084 0 : bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds
2085 :
2086 : ! The value for jj is negative to indicate that
2087 : ! the vector is from ii to jj
2088 0 : bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1
2089 0 : bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds
2090 :
2091 : ! The unitary vector is always from ii to jj
2092 0 : bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(rpt(:)-xcart(:,ii))/bl
2093 0 : bonds_tmp%bond_length(bonds_tmp%nbonds)=bl
2094 :
2095 0 : if (ii==jj) then
2096 0 : bonds_tmp%nbonds= bonds_tmp%nbonds+1
2097 : end if
2098 :
2099 : end if
2100 :
2101 : end do !! irshift
2102 :
2103 : end do !! jj
2104 : end do !! ii
2105 :
2106 6 : call print_bonds(amu,bonds_tmp,natom,ntypat,symbol,typat,znucl)
2107 :
2108 :
2109 : !write(std_out,*) 'make_bonds 05'
2110 : !##########################################################
2111 : !### 05. Deallocate all the arrays inside bonds
2112 : !### allocate them with the right size and fill them
2113 :
2114 6 : call bonds_free(bonds)
2115 :
2116 6 : bonds%nbonds=bonds_tmp%nbonds
2117 :
2118 6 : if (bonds%nbonds>0) then
2119 : ! Allocate the arrays with exactly the rigth nbonds
2120 18 : ABI_MALLOC(bonds%bond_vect,(3,bonds%nbonds))
2121 18 : ABI_MALLOC(bonds%bond_length,(bonds%nbonds))
2122 24 : ABI_MALLOC(bonds%indexi,(natom,bonds%nbonds))
2123 12 : ABI_MALLOC(bonds%nbondi,(natom))
2124 :
2125 : ! Fill the values
2126 30 : bonds%bond_vect(:,1:bonds%nbonds)=bonds_tmp%bond_vect(:,1:bonds%nbonds)
2127 12 : bonds%bond_length(1:bonds%nbonds)=bonds_tmp%bond_length(1:bonds%nbonds)
2128 24 : bonds%indexi(:,1:bonds%nbonds)=bonds_tmp%indexi(:,1:bonds%nbonds)
2129 18 : bonds%nbondi(:)=bonds_tmp%nbondi(:)
2130 : end if
2131 :
2132 6 : call bonds_free(bonds_tmp)
2133 :
2134 6 : end subroutine make_bonds_new
2135 : !!***
2136 :
2137 : !----------------------------------------------------------------------
2138 :
2139 : !!****f* m_abimover/bonds_free
2140 : !! NAME
2141 : !! bonds_free
2142 : !!
2143 : !! FUNCTION
2144 : !! Free memory
2145 : !!
2146 : !! SOURCE
2147 :
2148 18 : subroutine bonds_free(bonds)
2149 :
2150 : !Arguments ------------------------------------
2151 : type(go_bonds),intent(inout) :: bonds
2152 :
2153 : ! *********************************************************************
2154 :
2155 18 : ABI_SFREE(bonds%bond_vect)
2156 18 : ABI_SFREE(bonds%bond_length)
2157 18 : ABI_SFREE(bonds%nbondi)
2158 18 : ABI_SFREE(bonds%indexi)
2159 :
2160 18 : end subroutine bonds_free
2161 : !!***
2162 :
2163 : !----------------------------------------------------------------------
2164 :
2165 : !!****f* m_abimover/print_bonds
2166 : !! NAME
2167 : !! print_bonds
2168 : !!
2169 : !! FUNCTION
2170 : !! Print the bonds
2171 : !!
2172 : !! INPUTS
2173 : !! natom= Number of atoms
2174 : !! ntypat= Number of type of atoms
2175 : !! znucl= Z number of the atom
2176 : !!
2177 : !! OUTPUT
2178 : !! bonds= Structure that store all the information about
2179 : !! bonds created by this routine:
2180 : !! nbonds= Total number of bonds
2181 : !! bondij= Unitary vector along the bond direction
2182 : !! distij= Distances between atoms i and j (including shift)
2183 : !! listij= Indices of bonds going from i to j
2184 : !! listji= Indices of bonds going from j to i
2185 : !! indexij= Number of bonds between i and j
2186 : !! indexji= Number of bonds between j and i
2187 : !! tolerance
2188 : !!
2189 : !! SOURCE
2190 :
2191 6 : subroutine print_bonds(amu,bonds,natom,ntypat,symbol,typat,znucl)
2192 :
2193 : !Arguments ------------------------------------
2194 : !scalars
2195 : integer,intent(in) :: natom,ntypat
2196 : integer,intent(in) :: typat(natom)
2197 : real(dp),intent(in) :: znucl(ntypat)
2198 : real(dp),intent(in) :: amu(ntypat)
2199 : character(len=2),intent(in) :: symbol(ntypat)
2200 : type(go_bonds),intent(in) :: bonds
2201 :
2202 : !Local variables ------------------------------
2203 : !scalars
2204 : integer :: ii,jj,kk
2205 :
2206 : ! *********************************************************************
2207 :
2208 6 : write(std_out,'(a)') ch10
2209 438 : write(std_out,'(a,72a,a)') '---BONDS',('-',kk=1,72),ch10
2210 6 : write(std_out,'(a,i3)') 'Number of atoms: ',natom
2211 6 : write(std_out,'(a,i3)') 'Number of bonds: ',bonds%nbonds
2212 6 : write(std_out,'(a,f6.3,a,a)') 'Tolerance of bonds: ',bonds%tolerance,' times the sum of covalent radius',ch10
2213 :
2214 18 : do ii=1,natom
2215 12 : write(std_out,'(a,i3)') 'ATOM number: ',ii
2216 12 : write(std_out,'(a,f8.3)') ' Z: ',znucl(typat(ii))
2217 12 : write(std_out,'(a,f8.3)') ' Weight: ',amu(typat(ii))
2218 12 : write(std_out,'(a,a3)') ' Symbol: ',symbol(typat(ii))
2219 12 : write(std_out,'(a,i3)') ' Number of bonds: ',bonds%nbondi(ii)
2220 :
2221 30 : do jj=1,bonds%nbondi(ii)
2222 12 : write(std_out,'(a,i3,a,a,i3,a,3f7.3,a,f7.3)') ' [',jj,']',&
2223 12 : & ' Index of bond: ',bonds%indexi(ii,jj),&
2224 12 : & ' Unitary vector: ',bonds%bond_vect(:,abs(bonds%indexi(ii,jj))),&
2225 36 : & ' Bond length: ',bonds%bond_length(abs(bonds%indexi(ii,jj)))
2226 : end do
2227 :
2228 : end do
2229 :
2230 12 : do ii=1,bonds%nbonds
2231 6 : write(std_out,'(a,i3)') 'BOND Index=',ii
2232 6 : write(std_out,'(a,3f8.3)') ' Vector',bonds%bond_vect(:,ii)
2233 12 : write(std_out,'(a,f8.3)') ' bond Length',bonds%bond_length(ii)
2234 : end do
2235 :
2236 6 : end subroutine print_bonds
2237 : !!***
2238 :
2239 : !----------------------------------------------------------------------
2240 :
2241 : !!****f* m_abimover/delocint_ini
2242 : !!
2243 : !! NAME
2244 : !! delocint_ini
2245 : !!
2246 : !! FUNCTION
2247 : !! ini function for delocint object
2248 : !!
2249 : !! INPUT
2250 : !!
2251 : !! OUTPUT
2252 : !!
2253 : !! SIDE EFFECTS
2254 : !! deloc= container object for delocalized internal coordinates
2255 : !!
2256 : !! SOURCE
2257 :
2258 2 : subroutine delocint_ini(deloc)
2259 :
2260 : !Arguments ------------------------------------
2261 : !scalars
2262 : type(delocint), intent(out) :: deloc
2263 :
2264 : !Local variables ------------------------------
2265 : !scalars
2266 : integer :: ii,irshift,jj,kk,nshell
2267 :
2268 : ! *********************************************************************
2269 :
2270 2 : nshell=3
2271 2 : deloc%nrshift=(2*nshell+1)**3
2272 2 : deloc%icenter = nshell*(2*nshell+1)**2 + nshell*(2*nshell+1) + nshell + 1
2273 :
2274 2 : ABI_MALLOC(deloc%rshift,(3,deloc%nrshift))
2275 2 : irshift=0
2276 16 : do ii=-nshell,nshell
2277 114 : do jj=-nshell,nshell
2278 798 : do kk=-nshell,nshell
2279 686 : irshift=irshift+1
2280 2842 : deloc%rshift(:,irshift) = (/dble(ii),dble(jj),dble(kk)/)
2281 : end do
2282 : end do
2283 : end do
2284 :
2285 2 : end subroutine delocint_ini
2286 : !!***
2287 :
2288 : !----------------------------------------------------------------------
2289 :
2290 : !!****f* m_abimover/delocint_fin
2291 : !!
2292 : !! NAME
2293 : !! delocint_fin
2294 : !!
2295 : !! FUNCTION
2296 : !! destructor function for delocint object
2297 : !!
2298 : !! INPUT
2299 : !! deloc= container object for delocalized internal coordinates
2300 : !!
2301 : !! OUTPUT
2302 : !!
2303 : !! SOURCE
2304 :
2305 3976 : subroutine delocint_fin(deloc)
2306 :
2307 : type(delocint), intent(inout) :: deloc
2308 :
2309 3976 : ABI_SFREE(deloc%angs)
2310 3976 : ABI_SFREE(deloc%bonds)
2311 3976 : ABI_SFREE(deloc%carts)
2312 3976 : ABI_SFREE(deloc%dihedrals)
2313 3976 : ABI_SFREE(deloc%rshift)
2314 :
2315 3976 : end subroutine delocint_fin
2316 : !!***
2317 :
2318 : !----------------------------------------------------------------------
2319 :
2320 : !!****f* m_abimover/make_angles_new
2321 : !! NAME
2322 : !! make_angles_new
2323 : !!
2324 : !! FUNCTION
2325 : !! Fill the contents of the angles structure, that contains
2326 : !! all non redundant angles that could be generated between
2327 : !! all the atoms in the unitary cell and their adjacent cells
2328 : !! An angle is establish when an atom has two or more bonds.
2329 : !! The angles structure contains information about the atoms
2330 : !! involved, the value of the angle in radians, and the unitary
2331 : !! vector perpendicular to the plane of the three atoms that
2332 : !! build the angle.
2333 : !!
2334 : !! INPUTS
2335 : !! natom= Number of atoms
2336 : !! ntypat= Number of type of atoms
2337 : !! rprimd= Dimensional primitive vectors of the cell
2338 : !! xcart= Cartesian coordinates of the atoms
2339 : !! znucl= Z number of the atom
2340 : !! bonds= Structure that store all the information about
2341 : !! bonds created by this routine:
2342 : !! nbonds= Total number of bonds
2343 : !! nbondi= Number of bonds for atom i
2344 : !! indexi= Indeces of bonds for atom i
2345 : !! bond_length= Distances between atoms i and j (including shift)
2346 : !! bond_vect= Unitary vector for the bond from i to j
2347 : !! tolerance= The tolerance is multiplied to the
2348 : !! adition of covalent radius to decide if a bond is created
2349 : !!
2350 : !! OUTPUT
2351 : !! angles= Structure that store the information about
2352 : !! angles created by this routine
2353 : !! nangles= Total number of angles
2354 : !! angle_vertex= Index of the atom for that angle
2355 : !! angle_value= Value of the angle in radians
2356 : !! angle_bonds= Indices of the bonds
2357 : !! angle_vect= Unitary vector perpendicular to the plane of the angle
2358 : !!
2359 : !! SOURCE
2360 :
2361 : !This routine has been disables since it's broken
2362 : #if 0
2363 :
2364 : subroutine make_angles_new(angles,bonds,natom,ntypat,rprimd,typat,xcart,znucl)
2365 :
2366 : !Arguments ------------------------------------
2367 : !scalars
2368 : integer,intent(in) :: natom,ntypat
2369 : !arrays
2370 : integer,intent(in) :: typat(natom)
2371 : real(dp),intent(in) :: znucl(ntypat)
2372 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
2373 : type(go_bonds),intent(in) :: bonds
2374 : type(go_angles),intent(inout) :: angles
2375 :
2376 : !Local variables ------------------------------
2377 : !scalars
2378 : integer :: ii,jj,kk,iangle
2379 : type(atomdata_t) :: atom
2380 :
2381 : !arrays
2382 : type(go_bonds) :: bonds_tmp
2383 : character(len=2) :: symbol(ntypat)
2384 : real(dp) :: amu(ntypat)
2385 : integer :: shift(3,13) ! Represent all shift vectors that are not equivalent by central symmetry
2386 : ! For example (1,1,-1) is equivalent to (-1,-1,1)
2387 : ! It means that bond between atom i in the original cell and atom j in the
2388 : ! cell with cordinates (1,1,-1) is equivalent to the bond between atom j in
2389 : ! the orignal cell and atom i in the cell with coordinates (-1,-1,1)
2390 : ! The trivial shift (0,0,0) is excluded here
2391 : real(dp) :: rcov(ntypat) ! Covalent radius
2392 : real(dp) :: rpt(3)
2393 :
2394 : !***************************************************************************
2395 : !Beginning of executable session
2396 : !***************************************************************************
2397 : ABI_ERROR("This routine is not tested")
2398 :
2399 : !write(std_out,*) 'make_bonds 01'
2400 : !##########################################################
2401 : !### 01. Compute covalent radius
2402 :
2403 : do ii=1,ntypat
2404 : call atomdata_from_znucl(atom,znucl(ii))
2405 : amu(ii) = atom%amu
2406 : rcov(ii) = atom%rcov
2407 : symbol(ii) = symbol(ii)
2408 : end do
2409 :
2410 : !write(std_out,*) 'make_bonds 02'
2411 : !##########################################################
2412 : !### 02. Fill the 13 posible shift conecting adjacent cells
2413 :
2414 : shift(:,:)=reshape( (/ 1,0,0,&
2415 : & 0, 1, 0,&
2416 : & 0, 0, 1,&
2417 : & 1, 1, 0,&
2418 : & 1,-1, 0,&
2419 : & 0, 1, 1,&
2420 : & 0, 1,-1,&
2421 : & 1, 0, 1,&
2422 : & 1, 0,-1,&
2423 : & 1, 1, 1,&
2424 : & 1,-1, 1,&
2425 : & 1, 1,-1,&
2426 : & 1,-1,-1 /), (/ 3, 13 /))
2427 :
2428 : !write(std_out,*) 'make_bonds 03'
2429 : !##########################################################
2430 : !### 03. Initialize the values of bonds
2431 :
2432 : !The total number of bonds could not be predicted without
2433 : !compute all the distances, but the extreme case is linking
2434 : !all the atoms within all adjacent cells (natom*natom*13)
2435 : !plus the all the bonds inside the original cell (natom*(natom-1))
2436 :
2437 : bonds_tmp%nbonds=0
2438 : bonds_tmp%tolerance=bonds%tolerance
2439 : ibond=0
2440 :
2441 : ABI_MALLOC(bonds_tmp%bond_vect,(3,natom*natom*14-natom))
2442 : ABI_MALLOC(bonds_tmp%bond_length,(natom*natom*14-natom))
2443 :
2444 : !indexi contains the indices to the bonds
2445 : ABI_MALLOC(bonds_tmp%indexi,(natom,natom*natom*14-natom))
2446 :
2447 : ABI_MALLOC(bonds_tmp%nbondi,(natom))
2448 :
2449 : bonds_tmp%indexi(:,:)=0
2450 : bonds_tmp%nbondi(:)=0
2451 :
2452 : !write(std_out,*) 'make_bonds 04'
2453 : !##########################################################
2454 : !### 04. Compute the bonds inside the original cell
2455 : !### shift=(0,0,0)
2456 :
2457 : do ii=1,natom
2458 : rcov1 = rcov(typat(ii))
2459 :
2460 : do jj=ii+1,natom
2461 : rcov2 = rcov(typat(jj))
2462 :
2463 : bl=bond_length(xcart(:,ii),xcart(:,jj))
2464 :
2465 : if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then
2466 : ! We have a new bond, nbonds starts from
2467 : ! 0, so it could be used to index the
2468 : ! locations of bondij and distij
2469 :
2470 : ! Increase the number of bonds
2471 : bonds_tmp%nbonds= bonds_tmp%nbonds+1
2472 :
2473 : ! The number of bonds for atoms ii and jj
2474 : ! needs to raise by one
2475 : bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1
2476 : bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1
2477 :
2478 : bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds
2479 : ! The value for jj is negative to indicate that
2480 : ! the vector is from ii to jj
2481 : bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds
2482 :
2483 : ! The unitary vector is always from ii to jj
2484 : bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(xcart(:,jj)-xcart(:,ii))/bl
2485 : bonds_tmp%bond_length(bonds_tmp%nbonds)=bl
2486 :
2487 : end if
2488 :
2489 : end do !! jj
2490 : end do !! ii
2491 :
2492 : !write(std_out,*) 'make_bonds 05'
2493 : !##########################################################
2494 : !### 05. Compute the bonds outside the original cell
2495 : !### 13 shifts considered
2496 :
2497 : !Bonds between identical atoms but in diferent cells are
2498 : !allowed
2499 :
2500 : do ii=1,natom
2501 : rcov1 = rcov(typat(ii))
2502 : do jj=1,natom
2503 : rcov2 = rcov(typat(jj))
2504 :
2505 : do irshift=1,13
2506 :
2507 : do kk=1,3
2508 : rpt(kk) = xcart(kk,jj)+&
2509 : & shift(1,irshift)*rprimd(kk,1)+ &
2510 : & shift(2,irshift)*rprimd(kk,2)+ &
2511 : & shift(3,irshift)*rprimd(kk,3)
2512 : end do
2513 :
2514 :
2515 : bl =bond_length(xcart(:,ii),rpt)
2516 :
2517 : if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then
2518 :
2519 : ! We have a new bond, nbonds starts from
2520 : ! 0, so it could be used to index the
2521 : ! locations of bondij and distij
2522 :
2523 : ! Increase the number of bonds
2524 : bonds_tmp%nbonds= bonds_tmp%nbonds+1
2525 :
2526 : ! The number of bonds for atoms ii and jj
2527 : ! needs to raise by one
2528 : bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1
2529 : bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds
2530 :
2531 : ! The value for jj is negative to indicate that
2532 : ! the vector is from ii to jj
2533 : bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1
2534 : bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds
2535 :
2536 : ! The unitary vector is always from ii to jj
2537 : bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(rpt(:)-xcart(:,ii))/bl
2538 : bonds_tmp%bond_length(bonds_tmp%nbonds)=bl
2539 :
2540 : if (ii==jj) then
2541 : bonds_tmp%nbonds= bonds_tmp%nbonds+1
2542 : end if
2543 :
2544 : end if
2545 :
2546 : end do !! irshift
2547 :
2548 : end do !! jj
2549 : end do !! ii
2550 :
2551 : call print_bonds(amu,bonds_tmp,natom,ntypat,symbol,typat,znucl)
2552 :
2553 :
2554 : !write(std_out,*) 'make_bonds 05'
2555 : !##########################################################
2556 : !### 05. Deallocate all the arrays inside bonds
2557 : !### allocate them with the right size and fill them
2558 :
2559 : call bonds_free(bonds)
2560 :
2561 : bonds%nbonds=bonds_tmp%nbonds
2562 :
2563 : if (bonds%nbonds>0) then
2564 : ! Allocate the arrays with exactly the rigth nbonds
2565 : ABI_MALLOC(bonds%bond_vect,(3,bonds%nbonds))
2566 : ABI_MALLOC(bonds%bond_length,(bonds%nbonds))
2567 : ABI_MALLOC(bonds%indexi,(natom,bonds%nbonds))
2568 : ABI_MALLOC(bonds%nbondi,(natom))
2569 :
2570 : ! Fill the values
2571 : bonds%bond_vect(:,1:bonds%nbonds)=bonds_tmp%bond_vect(:,1:bonds%nbonds)
2572 : bonds%bond_length(1:bonds%nbonds)=bonds_tmp%bond_length(1:bonds%nbonds)
2573 : bonds%indexi(:,1:bonds%nbonds)=bonds_tmp%indexi(:,1:bonds%nbonds)
2574 : bonds%nbondi(:)=bonds_tmp%nbondi(:)
2575 : end if
2576 :
2577 : call bonds_free(bonds_tmp)
2578 :
2579 : end subroutine make_angles_new
2580 : !!***
2581 :
2582 : #endif
2583 :
2584 : !----------------------------------------------------------------------
2585 :
2586 0 : end module m_abimover
|
