Branch data 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 [ - + ]: 4403 : 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 : 4403 : natom=dtset%natom
445 : :
446 : 4403 : ab_mover%delayperm =dtset%delayperm
447 : 4403 : ab_mover%diismemory =dtset%diismemory
448 : 4403 : ab_mover%goprecon =dtset%goprecon
449 : 4403 : ab_mover%jellslab =dtset%jellslab
450 : 4403 : ab_mover%natom =dtset%natom
451 : 4403 : ab_mover%nconeq =dtset%nconeq
452 : 4403 : ab_mover%nnos =dtset%nnos
453 : 4403 : ab_mover%nsym =dtset%nsym
454 : 4403 : ab_mover%ntypat =dtset%ntypat
455 : 4403 : ab_mover%optcell =dtset%optcell
456 : 4403 : ab_mover%restartxf =dtset%restartxf
457 : 4403 : ab_mover%signperm =dtset%signperm
458 : 4403 : ab_mover%ionmov =dtset%ionmov
459 : 4403 : ab_mover%bmass =dtset%bmass
460 : 4403 : ab_mover%dtion =dtset%dtion
461 : 4403 : ab_mover%friction =dtset%friction
462 : 4403 : ab_mover%mdwall =dtset%mdwall
463 : 4403 : ab_mover%noseinert =dtset%noseinert
464 : 4403 : ab_mover%ph_nqshift =dtset%ph_nqshift
465 : 4403 : ab_mover%strprecon =dtset%strprecon
466 : 4403 : ab_mover%vis =dtset%vis
467 : :
468 : 4403 : ab_mover%iatfix =>dtset%iatfix(:,1:natom)
469 : 4403 : ab_mover%symafm =>dtset%symafm
470 : 4403 : ab_mover%symrel =>dtset%symrel
471 : 4403 : ab_mover%tnons =>dtset%tnons
472 : 4403 : ab_mover%ph_ngqpt =>dtset%ph_ngqpt
473 : 4403 : ab_mover%ph_qshift =>dtset%ph_qshift
474 : 4403 : ab_mover%typat =>dtset%typat(1:natom)
475 : 4403 : ab_mover%prtatlist =>dtset%prtatlist(1:natom)
476 : 4403 : ab_mover%goprecprm =>dtset%goprecprm
477 : 4403 : ab_mover%mdtemp =>dtset%mdtemp
478 : 4403 : ab_mover%strtarget =>dtset%strtarget
479 : 4403 : ab_mover%qmass =>dtset%qmass
480 : 4403 : ab_mover%znucl =>dtset%znucl
481 : :
482 : 4403 : ab_mover%use_pimd_routine = .false.
483 : :
484 : 4403 : ab_mover%amu_curr =>amu_curr
485 [ - + + - : 13209 : ABI_MALLOC(ab_mover%amass,(natom))
+ - + - -
+ - + ]
486 [ + - + + : 39504 : 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 [ + + + + : 65799 : ab_mover%ndof=count(ab_mover%iatfix==0)
+ + ]
492 : :
493 : : !Filename for Hessian matrix (NOT IN DTSET)
494 : 4403 : ab_mover%fnameabi_hes =>dtfil%fnameabi_hes
495 : : !Filename for _HIST file
496 : 4403 : 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 [ + + ]: 4403 : if (ab_mover%optcell/=0)then
509 : 215 : specs%isARused=.TRUE.
510 : : else
511 : 4188 : specs%isARused=.FALSE.
512 : : end if
513 : :
514 : : !Velocities are never changed excepts for ionmov=1,6,7,8
515 : 4403 : specs%isVused=.FALSE.
516 : :
517 : : !In general convergence is needed
518 : 4403 : specs%isFconv=.TRUE.
519 : :
520 : : !specs%ncycle is 1 by default except for ionmov=1,9,14
521 : 4403 : specs%ncycle=1
522 : :
523 : : !specs%nhist is -1 by default store all the history
524 : 4403 : 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 [ + + + + : 4403 : write(msg,"(a,i0)")"Wrong value for ionmov: ",ab_mover%ionmov
+ + + + +
+ + + + -
- + - + +
- - + ]
829 : : !ABI_ERROR(msg)
830 : : end select
831 : :
832 : 4403 : 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 : 4403 : subroutine abimover_destroy(ab_mover)
850 : :
851 : : !Arguments ------------------------------------
852 : : type(abimover),intent(inout) :: ab_mover
853 : :
854 : : ! ***************************************************************
855 : :
856 : 4403 : nullify(ab_mover%goprecprm)
857 : 4403 : nullify(ab_mover%iatfix)
858 : 4403 : nullify(ab_mover%mdtemp)
859 : 4403 : nullify(ab_mover%ph_ngqpt)
860 : 4403 : nullify(ab_mover%ph_qshift)
861 : :
862 : 4403 : nullify(ab_mover%prtatlist)
863 : 4403 : nullify(ab_mover%qmass)
864 : 4403 : nullify(ab_mover%strtarget)
865 : 4403 : nullify(ab_mover%symafm)
866 : 4403 : nullify(ab_mover%symrel)
867 : 4403 : nullify(ab_mover%tnons)
868 : 4403 : nullify(ab_mover%typat)
869 : 4403 : nullify(ab_mover%znucl)
870 : :
871 : 4403 : nullify(ab_mover%amu_curr)
872 [ - + ]: 4403 : ABI_FREE(ab_mover%amass)
873 : :
874 : 4403 : nullify(ab_mover%fnameabi_hes)
875 : 4403 : nullify(ab_mover%filnam_ds)
876 : :
877 : 4403 : 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 : 4041 : subroutine mttk_fin(mttk_vars)
1001 : :
1002 : : type(mttk_type), intent(inout) :: mttk_vars
1003 : :
1004 [ + + ]: 4041 : ABI_SFREE(mttk_vars%glogs)
1005 [ + + ]: 4041 : ABI_SFREE(mttk_vars%vlogs)
1006 [ + + ]: 4041 : ABI_SFREE(mttk_vars%xlogs)
1007 : :
1008 : 4041 : 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 : 4011 : subroutine delocint_fin(deloc)
2306 : :
2307 : : type(delocint), intent(inout) :: deloc
2308 : :
2309 [ + + ]: 4011 : ABI_SFREE(deloc%angs)
2310 [ + + ]: 4011 : ABI_SFREE(deloc%bonds)
2311 [ + + ]: 4011 : ABI_SFREE(deloc%carts)
2312 [ + + ]: 4011 : ABI_SFREE(deloc%dihedrals)
2313 [ + + ]: 4011 : ABI_SFREE(deloc%rshift)
2314 : :
2315 : 4011 : 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
# # # # #
# # # # #
# # # # #
# # # ]
|
