TABLE OF CONTENTS
- ABINIT/m_abimover
- defs_mover/abimover_print
- m_abimover/ab_xfh_type
- m_abimover/abiforstr
- m_abimover/abiforstr_fin
- m_abimover/abiforstr_ini
- m_abimover/abimover
- m_abimover/abimover_destroy
- m_abimover/abimover_ini
- m_abimover/abimover_specs
- m_abimover/angle_ang
- m_abimover/angle_dihedral
- m_abimover/bond_length
- m_abimover/bonds_free
- m_abimover/calc_prim_int
- m_abimover/delocint
- m_abimover/delocint_fin
- m_abimover/delocint_ini
- m_abimover/go_angles
- m_abimover/go_bonds
- m_abimover/make_angles
- m_abimover/make_angles_new
- m_abimover/make_bonds
- m_abimover/make_bonds_new
- m_abimover/make_dihedrals
- m_abimover/make_prim_internals
- m_abimover/mttk_fin
- m_abimover/mttk_ini
- m_abimover/mttk_type
- m_abimover/print_bonds
ABINIT/m_abimover [ Modules ]
NAME
m_abimover
FUNCTION
This module contains definition the types abimover, mttk, abiforstr, delocint, and bonds and their related ini and free routines
COPYRIGHT
Copyright (C) 2001-2024 ABINIT group (DCA, XG, GMR, SE, Mver, JJ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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
defs_mover/abimover_print [ Functions ]
NAME
abimover_print
FUNCTION
Print all the variables in a ab_mover
INPUTS
OUTPUT
SIDE EFFECTS
ab_mover <type(abimover)> = The ab_mover to nullify
NOTES
At present 29 variables are present in abimover if a new variable is added in abimover it should be added also for print here
SOURCE
867 subroutine abimover_print(ab_mover,iout) 868 869 !Arguments ------------------------------------ 870 integer,intent(in) :: iout 871 type(abimover),intent(inout) :: ab_mover 872 873 !Local variables------------------------------- 874 !arrays 875 character(len=1200) :: message 876 character(len=110) :: fmt 877 878 ! *************************************************************** 879 880 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)' 881 882 write(message,fmt)& 883 & 'Delta Time for IONs',ab_mover%dtion,ch10, & 884 & 'include a JELLium SLAB in the cell',ab_mover%jellslab,ch10, & 885 & 'Number of ATOMs',ab_mover%natom,ch10, & 886 & 'Number of CONstraint EQuations',ab_mover%nconeq,ch10, & 887 & 'Number of SYMmetry operations',ab_mover%nsym,ch10, & 888 & 'OPTimize the CELL shape and dimensions',ab_mover%optcell,ch10, & 889 & 'RESTART Xcart and Fred',ab_mover%restartxf,ch10, & 890 & 'Molecular Dynamics Initial Temperature',ab_mover%mdtemp(1),ch10, & 891 & 'Molecular Dynamics Final Temperature',ab_mover%mdtemp(2),ch10, & 892 & 'NOSE thermostat INERTia factor',ab_mover%noseinert,ch10, & 893 & 'STRess PRECONditioner',ab_mover%strprecon,ch10, & 894 & 'VIScosity',ab_mover%vis,ch10 895 896 ! ! arrays 897 ! ! Indices of AToms that are FIXed 898 ! integer, pointer :: iatfix(:,:) 899 ! ! SYMmetries, Anti-FerroMagnetic characteristics 900 ! integer, pointer :: symafm(:) 901 ! ! SYMmetry in REaL space 902 ! integer, pointer :: symrel(:,:,:) 903 ! Translation NON-Symmorphic vectors 904 ! real(dp), pointer :: tnons(:,:) 905 ! ! Mass of each atom (NOT IN DTSET) 906 ! real(dp), pointer :: amass(:) 907 ! ! STRess TARGET 908 ! real(dp), pointer :: strtarget(:) 909 ! Filename for Hessian matrix 910 ! character(len=fnlen), pointer :: fnameabi_hes 911 912 write(iout,*) 'CONTENT of ab_mover (scalar only)' 913 write(iout,'(a)') message 914 915 end subroutine abimover_print
m_abimover/ab_xfh_type [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
ab_xfh_type
FUNCTION
Datatype with the old structure for storing history used in gstate and brdmin,delocint, and others
NOTES
This is a transitional structure, to bridge between the old code and the new one base on abihist
SOURCE
318 type, public :: ab_xfh_type 319 320 integer :: nxfh,nxfhr,mxfh 321 ! mxfh = last dimension of the xfhist array 322 ! nxfh = actual number of (x,f) history pairs, see xfhist array 323 324 real(dp),allocatable :: xfhist(:,:,:,:) 325 ! xfhist(3,natom+4,2,mxfh) = (x,f) history array, also including rprim and stress 326 327 end type ab_xfh_type
m_abimover/abiforstr [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
abiforstr
FUNCTION
Store forces, stress and energy, cartesian and reduced forces one scalar for energy and 6 element array for stress
NOTES
SOURCE
285 type, public :: abiforstr 286 287 ! scalars 288 real(dp) :: etotal 289 ! Total energy 290 291 ! arrays 292 real(dp),allocatable :: fcart(:,:) 293 ! Cartesian forces 294 real(dp),allocatable :: gred(:,:) 295 ! Reduced forces 296 real(dp) :: strten(6) 297 ! Stress tensor (Symmetrical 3x3 matrix) 298 299 end type abiforstr
m_abimover/abiforstr_fin [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
abiforstr_fin
FUNCTION
destructor function for abiforstr object INPUT forstr
OUTPUT
SOURCE
1020 subroutine abiforstr_fin(forstr) 1021 1022 type(abiforstr), intent(inout) :: forstr 1023 1024 ABI_SFREE(forstr%fcart) 1025 ABI_SFREE(forstr%gred) 1026 1027 end subroutine abiforstr_fin
m_abimover/abiforstr_ini [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
abiforstr_ini
FUNCTION
destructor function for abiforstr object INPUT forstr
OUTPUT
SOURCE
992 subroutine abiforstr_ini(forstr,natom) 993 994 integer,intent(in) :: natom 995 type(abiforstr), intent(out) :: forstr 996 997 ABI_MALLOC(forstr%fcart,(3,natom)) 998 ABI_MALLOC(forstr%gred,(3,natom)) 999 1000 end subroutine abiforstr_ini
m_abimover/abimover [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
abimover
FUNCTION
This datatype has the purpose of store all the data taken usually from dtset (but not only) needed for the different predictors to update positions, acell, etc.
SOURCE
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 Shifts for the Qpoint Grid (used for ionmov 26 and 27) 87 integer :: ph_nqshift 88 ! Use by pred_isothermal only 89 integer :: nnos 90 ! Number of SYMmetry operations 91 integer :: nsym 92 ! Number of Types of atoms 93 integer :: ntypat 94 ! OPTimize the CELL shape and dimensions 95 integer :: optcell 96 ! RESTART Xcart and Fred 97 integer :: restartxf 98 ! Sign of Permutation (Used by pred_langevin only) 99 integer :: signperm 100 ! Ion movement 101 integer :: ionmov 102 103 ! Use by pred_isothermal only 104 real(dp) :: bmass 105 ! Delta Time for IONs 106 real(dp) :: dtion 107 ! Used by pred_langevin only 108 real(dp) :: friction 109 ! Used by pred_langevin only 110 real(dp) :: mdwall 111 ! Used by pred_nose only 112 real(dp) :: noseinert 113 ! STRess PRECONditioner 114 real(dp) :: strprecon 115 ! VIScosity 116 real(dp) :: vis 117 118 ! arrays 119 ! Indices of AToms that are FIXed 120 integer,pointer :: iatfix(:,:) ! iatfix(3,natom) 121 ! SYMmetries, Anti-FerroMagnetic characteristics 122 integer,pointer :: symafm(:) ! symafm(nsym) 123 ! SYMmetry in REaL space 124 integer,pointer :: symrel(:,:,:) ! symrel(3,3,nsym) 125 ! Translation NON-Symmorphic vectors 126 real(dp),pointer :: tnons(:,:) ! tnons(3,nsym) 127 ! TYPe of ATom 128 integer,pointer :: typat(:) ! typat(natom) 129 ! PRTint ATom LIST 130 integer,pointer :: prtatlist(:) ! prtatlist(natom) 131 ! Qpoint grid (used for ionmov 26 and 27) 132 integer,pointer :: ph_ngqpt(:) ! ph_ngqpt(3) 133 ! shift of the Qpoint Grid (used for ionmov 26 and 27) 134 real(dp),pointer :: ph_qshift(:,:) ! 135 ! amu input var for the current image 136 real(dp), pointer :: amu_curr(:) ! amu_curr(ntypat) 137 ! Mass of each atom 138 real(dp),pointer :: amass(:) ! amass(natom) 139 ! Geometry Optimization Preconditioner PaRaMeters 140 real(dp),pointer :: goprecprm(:) 141 ! Molecular Dynamics Initial and Final Temperature 142 real(dp),pointer :: mdtemp(:) ! mdtemp(2) (initial,final) 143 ! STRess TARGET 144 real(dp),pointer :: strtarget(:) ! strtarget(6) 145 ! Use by pred_isothermal only 146 real(dp),pointer :: qmass(:) 147 ! Z number of each NUCLeus 148 real(dp),pointer :: znucl(:) ! znucl(npsp) 149 150 ! Filename for Hessian matrix 151 character(len=fnlen), pointer :: fnameabi_hes 152 ! Filename for _HIST file 153 character(len=fnlen), pointer :: filnam_ds(:) ! dtfil%filnam_ds(5) 154 155 end type abimover
m_abimover/abimover_destroy [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
abimover_destroy
FUNCTION
Destroy the abimover structure
SIDE EFFECTS
ab_mover <type(abimover)> = The abimover structure to be destroyed
SOURCE
813 subroutine abimover_destroy(ab_mover) 814 815 !Arguments ------------------------------------ 816 type(abimover),intent(inout) :: ab_mover 817 818 ! *************************************************************** 819 820 nullify(ab_mover%goprecprm) 821 nullify(ab_mover%iatfix) 822 nullify(ab_mover%mdtemp) 823 nullify(ab_mover%ph_ngqpt) 824 nullify(ab_mover%ph_qshift) 825 826 nullify(ab_mover%prtatlist) 827 nullify(ab_mover%qmass) 828 nullify(ab_mover%strtarget) 829 nullify(ab_mover%symafm) 830 nullify(ab_mover%symrel) 831 nullify(ab_mover%tnons) 832 nullify(ab_mover%typat) 833 nullify(ab_mover%znucl) 834 835 nullify(ab_mover%amu_curr) 836 ABI_FREE(ab_mover%amass) 837 838 nullify(ab_mover%fnameabi_hes) 839 nullify(ab_mover%filnam_ds) 840 841 end subroutine abimover_destroy
m_abimover/abimover_ini [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
abimover_ini
FUNCTION
Initializes the abimover structure and the abimover_specs information
INPUTS
OUTPUT
NOTES
SOURCE
416 subroutine abimover_ini(ab_mover,amu_curr,dtfil,dtset,specs) 417 418 !Arguments ------------------------------------ 419 real(dp),target, intent(in) :: amu_curr(:) ! amu_curr(ntype) 420 type(abimover),intent(out) :: ab_mover 421 type(datafiles_type),target,intent(in) :: dtfil 422 type(dataset_type),target,intent(in) :: dtset 423 type(abimover_specs),intent(out) :: specs 424 425 !Local variables------------------------------- 426 !scalars 427 integer :: iatom,natom 428 character(len=500) :: msg 429 !arrays 430 431 ! *************************************************************** 432 433 434 !write(std_out,*) 'mover 01' 435 !########################################################### 436 !### 01. Initialization of ab_mover 437 438 !Copy or create pointers for the information from the Dataset (dtset) to the ab_mover structure 439 natom=dtset%natom 440 441 ab_mover%delayperm =dtset%delayperm 442 ab_mover%diismemory =dtset%diismemory 443 ab_mover%goprecon =dtset%goprecon 444 ab_mover%jellslab =dtset%jellslab 445 ab_mover%natom =dtset%natom 446 ab_mover%nconeq =dtset%nconeq 447 ab_mover%nnos =dtset%nnos 448 ab_mover%nsym =dtset%nsym 449 ab_mover%ntypat =dtset%ntypat 450 ab_mover%optcell =dtset%optcell 451 ab_mover%restartxf =dtset%restartxf 452 ab_mover%signperm =dtset%signperm 453 ab_mover%ionmov =dtset%ionmov 454 ab_mover%bmass =dtset%bmass 455 ab_mover%dtion =dtset%dtion 456 ab_mover%friction =dtset%friction 457 ab_mover%mdwall =dtset%mdwall 458 ab_mover%noseinert =dtset%noseinert 459 ab_mover%ph_nqshift =dtset%ph_nqshift 460 ab_mover%strprecon =dtset%strprecon 461 ab_mover%vis =dtset%vis 462 463 ab_mover%iatfix =>dtset%iatfix(:,1:natom) 464 ab_mover%symafm =>dtset%symafm 465 ab_mover%symrel =>dtset%symrel 466 ab_mover%tnons =>dtset%tnons 467 ab_mover%ph_ngqpt =>dtset%ph_ngqpt 468 ab_mover%ph_qshift =>dtset%ph_qshift 469 ab_mover%typat =>dtset%typat(1:natom) 470 ab_mover%prtatlist =>dtset%prtatlist(1:natom) 471 ab_mover%goprecprm =>dtset%goprecprm 472 ab_mover%mdtemp =>dtset%mdtemp 473 ab_mover%strtarget =>dtset%strtarget 474 ab_mover%qmass =>dtset%qmass 475 ab_mover%znucl =>dtset%znucl 476 477 ab_mover%amu_curr =>amu_curr 478 ABI_MALLOC(ab_mover%amass,(natom)) 479 do iatom=1,natom 480 ab_mover%amass(iatom)=amu_emass*amu_curr(dtset%typat(iatom)) 481 end do 482 483 !Filename for Hessian matrix (NOT IN DTSET) 484 ab_mover%fnameabi_hes =>dtfil%fnameabi_hes 485 !Filename for _HIST file 486 ab_mover%filnam_ds =>dtfil%filnam_ds 487 488 !call abimover_print(ab_mover,ab_out) 489 490 !write(std_out,*) 'mover 02' 491 !########################################################### 492 !### 02. Particularities of each predictor 493 494 !Default values first 495 !-------------------- 496 497 !acell and rprimd are never changed except if optcell/=0 498 if (ab_mover%optcell/=0)then 499 specs%isARused=.TRUE. 500 else 501 specs%isARused=.FALSE. 502 end if 503 504 !Velocities are never changed excepts for ionmov=1,6,7,8 505 specs%isVused=.FALSE. 506 507 !In general convergence is needed 508 specs%isFconv=.TRUE. 509 510 !specs%ncycle is 1 by default except for ionmov=1,9,14 511 specs%ncycle=1 512 513 !specs%nhist is -1 by default store all the history 514 specs%nhist=-1 515 516 !This is the initialization for ionmov==1 517 !----------------------------------------- 518 select case (ab_mover%ionmov) 519 case (1) 520 specs%ncycle=4 ! Number of internal cycles for first itime 521 specs%isFconv=.FALSE. ! Convergence is not used for MD 522 specs%isVused=.TRUE. ! Velocities are used 523 ! TEMPORARLY optcell is not allow 524 specs%isARused=.FALSE. 525 ! Values use in XML Output 526 specs%type4xml='moldyn' 527 specs%crit4xml='none' 528 ! Name of specs%method 529 if (abs(ab_mover%vis)<=1.d-8) then 530 specs%method = 'Molecular dynamics without viscosity (vis=0)' 531 else 532 write(specs%method,'(a,1p,e12.5,a)')& 533 'Molecular dynamics with viscosity (vis=',ab_mover%vis,')' 534 end if 535 ! Number of history 536 specs%nhist = 6 537 ! This is the initialization for ionmov==2,3 538 ! ------------------------------------------- 539 case (2,3) 540 ! Values use in XML Output 541 specs%type4xml='bfgs' 542 specs%crit4xml='tolmxf' 543 ! Name of specs%method 544 if (ab_mover%ionmov==2) then 545 specs%method = 'Broyden-Fletcher-Goldfard-Shanno method (forces)' 546 else 547 specs%method = 'Broyden-Fletcher-Goldfard-Shanno method (forces,Tot energy)' 548 end if 549 ! Number of history 550 specs%nhist = 3 551 ! This is the initialization for ionmov==4,5 552 ! ------------------------------------------- 553 case (4,5) 554 ! Values used in XML Output 555 specs%type4xml='simple' 556 specs%crit4xml='tolmxf' 557 ! Name of specs%method 558 if (ab_mover%ionmov==4) then 559 specs%method = 'Conjugate gradient of potential and ionic degrees of freedom' 560 else 561 specs%method = 'Simple relaxation of ionic positions' 562 end if 563 ! Number of history 564 specs%nhist = 3 565 ! This is the initialization for ionmov==6 566 ! ------------------------------------------ 567 case (6) 568 specs%isFconv=.FALSE. ! Convergence is not used for MD 569 ! TEMPORARLY optcell is not allow 570 specs%isARused=.FALSE. 571 specs%isVused=.TRUE. ! Velocities are used 572 ! Values use in XML Output 573 specs%type4xml='verlet' 574 specs%crit4xml='tolmxf' 575 ! Name of specs%method 576 specs%method = 'Verlet algorithm for molecular dynamics' 577 ! Number of history 578 specs%nhist = 3 579 ! This is the initialization for ionmov==7 580 ! ------------------------------------------ 581 case (7) 582 ! TEMPORARLY optcell is not allow 583 specs%isARused=.FALSE. 584 specs%isVused=.TRUE. ! Velocities are used 585 ! Values use in XML Output 586 specs%type4xml='verlet' 587 specs%crit4xml='tolmxf' 588 ! Name of specs%method 589 specs%method = 'Verlet algorithm blocking every atom where dot(vel,force)<0' 590 ! Number of history 591 specs%nhist = 3 592 ! This is the initialization for ionmov==8 593 ! ------------------------------------------ 594 case (8) 595 specs%isVused=.TRUE. 596 ! TEMPORARLY optcell is not allow 597 specs%isARused=.FALSE. 598 ! Values use in XML Output 599 specs%type4xml='nose' 600 specs%crit4xml='tolmxf' 601 ! Name of specs%method 602 specs%method = 'Verlet algorithm with a nose-hoover thermostat' 603 ! Number of history 604 specs%nhist = 3 605 ! This is the initialization for ionmov==9 606 ! ------------------------------------------ 607 case (9) 608 ! TEMPORARLY optcell is not allow 609 specs%isARused=.FALSE. 610 specs%isVused=.TRUE. ! Velocities are used 611 specs%ncycle=3 612 ! Values use in XML Output 613 specs%type4xml='langevin' 614 specs%crit4xml='tolmxf' 615 ! Name of specs%method 616 specs%method = 'Langevin molecular dynamics' 617 ! Number of history 618 specs%nhist = 3 619 ! This is the initialization for ionmov==10 and 11 620 ! ------------------------------------------- 621 case (10,11) 622 ! TEMPORARLY optcell is not allow 623 specs%isARused=.FALSE. 624 ! Values use in XML Output 625 if(ab_mover%ionmov==10)specs%type4xml='delocint' 626 if(ab_mover%ionmov==11)specs%type4xml='cg' 627 specs%crit4xml='tolmxf' 628 ! Name of specs%method 629 if(ab_mover%ionmov==10)specs%method = 'BFGS with delocalized internal coordinates' 630 if(ab_mover%ionmov==11)specs%method = 'Conjugate gradient with deloc. int. coord.' 631 ! Number of history 632 specs%nhist = 3 633 ! This is the initialization for ionmov==12 634 ! ------------------------------------------- 635 case (12) 636 ! TEMPORARLY optcell is not allow 637 specs%isARused=.FALSE. 638 specs%isVused=.TRUE. ! Velocities are used 639 ! Values use in XML Output 640 specs%isFconv=.FALSE. ! Convergence is not used for MD 641 specs%type4xml='isokin' 642 specs%crit4xml='tolmxf' 643 ! Name of specs%method 644 specs%method = 'Isokinetic ensemble molecular dynamics' 645 ! Number of history 646 specs%nhist = 3 647 ! This is the initialization for ionmov==13 648 ! ------------------------------------------- 649 case (13) 650 ! optcell is allow 651 specs%isARused=.TRUE. ! RPRIMD and ACELL may change 652 specs%isVused=.TRUE. ! Velocities are used 653 specs%isFconv=.FALSE. ! Convergence is not used for MD 654 ! Values use in XML Output 655 specs%type4xml='isother' 656 specs%crit4xml='tolmxf' 657 ! Name of specs%method 658 specs%method = 'Isothermal/isenthalpic ensemble molecular dynamics' 659 ! Number of history 660 specs%nhist = 3 661 ! This is the initialization for ionmov==14 662 ! ------------------------------------------- 663 case (14) 664 specs%ncycle=16 665 specs%isFconv=.FALSE. ! Convergence is not used for MD 666 specs%isVused=.TRUE. ! Velocities are used 667 ! TEMPORARLY optcell is not allow 668 specs%isARused=.FALSE. 669 ! Values use in XML Output 670 specs%type4xml='srkna14' 671 specs%crit4xml='tolmxf' 672 ! Name of specs%method 673 specs%method = 'Symplectic algorithm Runge-Kutta-Nystrom SRKNa14' 674 ! Number of history 675 specs%nhist = 3 676 677 ! This is the initialization for ionmov==15 678 ! ------------------------------------------- 679 case (15) 680 ! Values use in XML Output 681 specs%type4xml='FIRE' 682 specs%isVused=.TRUE. ! Velocities are used 683 specs%isARused=.TRUE. 684 specs%crit4xml='tolmxf' 685 ! Name of specs%method 686 specs%method = 'Fast inertial relaxation engine' 687 ! Number of history 688 specs%nhist = 2 689 ! This is the initialization for ionmov==20 690 ! ------------------------------------------- 691 case (20) 692 ! TEMPORARLY optcell is not allow 693 specs%isARused=.FALSE. 694 ! Values use in XML Output 695 specs%type4xml='diisrelax' 696 specs%crit4xml='tolmxf' 697 ! Name of specs%method 698 specs%method = 'Ionic positions relaxation using DIIS' 699 ! Number of history 700 specs%nhist = 3 701 ! This is the initialization for ionmov==21 702 ! ------------------------------------------- 703 case (21) 704 specs%isARused=.TRUE. 705 ! Values use in XML Output 706 specs%type4xml='steepdesc' 707 specs%crit4xml='tolmxf' 708 ! Name of specs%method 709 specs%method = 'Steepest descend algorithm' 710 ! Number of history 711 specs%nhist = 3 712 ! This is the initialization for ionmov==22 713 ! ------------------------------------------- 714 case (22) 715 ! Values use in XML Output 716 specs%type4xml='lbfgs' 717 specs%crit4xml='tolmxf' 718 ! Name of specs%method 719 specs%method = 'Limited-memory Broyden-Fletcher-Goldfard-Shanno method' 720 ! Number of history 721 specs%nhist = 3 722 ! This is the initialization for ionmov==23 723 ! ------------------------------------------- 724 case (23) 725 specs%ncycle=2 726 ! TEMPORARLY optcell is not allow 727 specs%isARused=.FALSE. 728 specs%isVused=.TRUE. ! Velocities are used 729 ! Values use in XML Output 730 specs%isFconv=.FALSE. ! Convergence is not used for MD 731 specs%type4xml='isokin' 732 specs%crit4xml='tolmxf' 733 ! Name of specs%method 734 specs%method = 'Using LOTF Molecular dynamics' 735 ! Number of history 736 specs%nhist = 3 737 ! This is the initialization for ionmov==24 738 ! ------------------------------------------- 739 case (24) 740 specs%ncycle=1 741 ! TEMPORARLY optcell is not allow 742 specs%isARused=.FALSE. 743 specs%isVused=.TRUE. ! Velocities are used 744 ! Values use in XML Output 745 specs%isFconv=.FALSE. ! Convergence is not used for MD 746 specs%type4xml='velver' 747 specs%crit4xml='none' 748 ! Name of specs%method 749 specs%method = 'Symplectic velocity verlet Molecular dynamics' 750 ! Number of history 751 specs%nhist = 3 752 ! This is the initialization for ionmov==25 753 ! ------------------------------------------- 754 case (25) ! Hybrid Monte Carlo algorithm (fixed lattice vectors) 755 specs%ncycle = 12 ! Number of internal cycles (10+2) 756 specs%isFconv=.FALSE. ! Convergence is not used for Monte Carlo 757 specs%isVused=.TRUE. ! Velocities are used for update of atomic positions 758 ! optcell is not allowed 759 specs%isARused=.FALSE. 760 ! Values use in XML Output 761 specs%type4xml='hmc' 762 specs%crit4xml='none' 763 ! Name of specs%method 764 specs%method = 'Hybrid Monte Carlo' 765 specs%nhist = 3 766 ! This is the initialization for ionmov==27 767 ! ------------------------------------------- 768 case (27) ! Generation of the training set for effective potential 769 specs%ncycle = 1 ! Number of internal cycles 770 specs%isFconv=.FALSE. ! Convergence is not used 771 specs%isVused=.FALSE. ! Velocities are not used for update of atomic positions 772 ! Values use in XML Output 773 specs%type4xml='TS' 774 specs%crit4xml='none' 775 ! Name of specs%method 776 specs%method = 'training set generator' 777 ! Number of history 778 specs%nhist = -1 779 780 case (28) 781 ! Values used in XML Output 782 specs%type4xml='simple' 783 specs%crit4xml='tolmxf' 784 ! Name of specs%method 785 specs%method = "i-pi protocol" 786 ! Number of history 787 !specs%nhist = 3 788 ! This is the initialization for ionmov==6 789 790 case default 791 ! MG TODO: Why is this check deactivated. We should have an empty case for all the ionmov that can use the default values 792 write(msg,"(a,i0)")"Wrong value for ionmov: ",ab_mover%ionmov 793 !ABI_ERROR(msg) 794 end select 795 796 end subroutine abimover_ini
m_abimover/abimover_specs [ Types ]
[ Top ] [ m_abimover ] [ Types ]
m_abimover/angle_ang [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
angle_ang
FUNCTION
INPUTS
OUTPUT
SOURCE
1746 pure function angle_ang(r1,r2,r3) 1747 1748 !Arguments ------------------------------------ 1749 !scalars 1750 real(dp) :: angle_ang 1751 !arrays 1752 real(dp),intent(in) :: r1(3),r2(3),r3(3) 1753 1754 !Local variables ------------------------------ 1755 !scalars 1756 real(dp) :: cos_ang,n1,n2 1757 !arrays 1758 real(dp) :: rpt12(3),rpt32(3) 1759 1760 !****************************************************************** 1761 n1=bond_length(r1,r2) 1762 n2=bond_length(r3,r2) 1763 1764 rpt12(:) = r1(:)-r2(:) 1765 rpt32(:) = r3(:)-r2(:) 1766 1767 cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2 1768 1769 if (cos_ang > one - epsilon(one)*two) then 1770 cos_ang = one 1771 else if(cos_ang < -one + epsilon(one)*two) then 1772 cos_ang = -one 1773 end if 1774 1775 angle_ang=acos(cos_ang) 1776 1777 end function angle_ang
m_abimover/angle_dihedral [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
angle_dihedral
FUNCTION
INPUTS
OUTPUT
SOURCE
1792 function angle_dihedral(r1,r2,r3,r4) 1793 1794 !Arguments ------------------------------------ 1795 !scalars 1796 real(dp) :: angle_dihedral 1797 !arrays 1798 real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3) 1799 1800 !Local variables------------------------------------ 1801 !scalars 1802 real(dp) :: cos_dihedral,dih_sign,n1,n2,sin_dihedral 1803 !arrays 1804 real(dp) :: cp1232(3),cp3432(3),cpcp(3),rpt12(3),rpt32(3),rpt34(3) 1805 1806 !****************************************************************** 1807 1808 rpt12(:) = r1(:)-r2(:) 1809 rpt32(:) = r3(:)-r2(:) 1810 rpt34(:) = r3(:)-r4(:) 1811 1812 call acrossb(rpt12,rpt32,cp1232) 1813 call acrossb(rpt34,rpt32,cp3432) 1814 1815 !DEBUG 1816 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232 1817 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432 1818 !ENDDEBUG 1819 1820 n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2) 1821 n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2) 1822 1823 cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2 1824 !we use complementary of standard angle, so 1825 cos_dihedral = -cos_dihedral 1826 1827 call acrossb(cp1232,cp3432,cpcp) 1828 cpcp(:) = cpcp(:)/n1/n2 1829 sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))& 1830 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2) 1831 dih_sign = one 1832 !if (abs(sin_dihedral) > tol12) then 1833 !dih_sign = sin_dihedral/abs(sin_dihedral) 1834 !end if 1835 if (sin_dihedral < -tol12) then 1836 dih_sign = -one 1837 end if 1838 1839 !DEBUG 1840 !write(std_out,'(a,3E20.10)') 'angle_dihedral : cos sin dih_sign= ',& 1841 !& cos_dihedral,sin_dihedral,dih_sign 1842 !ENDDEBUG 1843 1844 if (cos_dihedral > one - epsilon(one)*two) then 1845 cos_dihedral = one 1846 else if(cos_dihedral < -one + epsilon(one)*two) then 1847 cos_dihedral = -one 1848 end if 1849 1850 angle_dihedral = dih_sign*acos(cos_dihedral) 1851 1852 end function angle_dihedral
m_abimover/bond_length [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
bond_length
FUNCTION
INPUTS
OUTPUT
SOURCE
1715 pure function bond_length(r1,r2) 1716 1717 !Arguments ------------------------------------ 1718 !scalars 1719 real(dp) :: bond_length 1720 !arrays 1721 real(dp),intent(in) :: r1(3),r2(3) 1722 1723 !Local variables ------------------------------------ 1724 !arrays 1725 real(dp) :: rpt(3) 1726 1727 !****************************************************************** 1728 rpt(:) = r1(:)-r2(:) 1729 bond_length = sqrt(rpt(1)**2+rpt(2)**2+rpt(3)**2) 1730 1731 end function bond_length
m_abimover/bonds_free [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
bonds_free
FUNCTION
Free memory
SOURCE
2117 subroutine bonds_free(bonds) 2118 2119 !Arguments ------------------------------------ 2120 type(go_bonds),intent(inout) :: bonds 2121 2122 ! ********************************************************************* 2123 2124 ABI_SFREE(bonds%bond_vect) 2125 ABI_SFREE(bonds%bond_length) 2126 ABI_SFREE(bonds%nbondi) 2127 ABI_SFREE(bonds%indexi) 2128 2129 end subroutine bonds_free
m_abimover/calc_prim_int [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
calc_prim_int
FUNCTION
calculate values of primitive internal coordinates as a function of cartesian ones.
INPUTS
angs= number of angles
bonds(2,2,nbond)=for a bond between iatom and jatom
bonds(1,1,nbond) = iatom
bonds(2,1,nbond) = icenter
bonds(1,2,nbond) = jatom
bonds(2,2,nbond) = irshift
carts(2,ncart)= index of total primitive internal, and atom (carts(2,:))
dihedrals(2,4,ndihed)=indexes to characterize dihedrals
dtset <type(dataset_type)>=all input variables for this dataset
nang(2,3,nang)=indexes to characterize angles
nbond=number of bonds
ncart=number of cartesian coordinates used
ndihed= number of dihedrals
ninternal=nbond+nang+ndihed+ncart: number of internal coordinates
nrshift= dimension of rshift
rprimd(3,3)=dimensional real space primitive translations (bohr)
rshift(3,nrshift)=shift in xred that must be done to find all neighbors of
a given atom within a given number of neighboring shells
xcart(3,natom)=cartesian coordinates of atoms (bohr)
OUTPUT
prim_int(ninternal)=values of primitive internal coordinates
SIDE EFFECTS
NOTES
SOURCE
1590 subroutine calc_prim_int(deloc,natom,rprimd,xcart,prim_int) 1591 1592 !Arguments ------------------------------------ 1593 !scalars 1594 integer,intent(in) :: natom 1595 type(delocint),intent(in) :: deloc 1596 !arrays 1597 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 1598 real(dp),intent(out) :: prim_int(deloc%ninternal) 1599 1600 !Local variables ------------------------------ 1601 !scalars 1602 integer :: i1,i2,i3,i4,iang,ibond,icart,idihed,iprim,s1,s2,s3,s4 1603 !arrays 1604 real(dp) :: r1(3),r2(3),r3(3),r4(3) 1605 1606 !************************************************************************ 1607 1608 !DEBUG 1609 !write(std_out,*) ' calc_prim_int : enter' 1610 !write(std_out,*) shape(deloc%bonds) 1611 !do ibond=1,deloc%nbond 1612 !do i1=1,2 1613 !write(std_out,'(2I5)') deloc%bonds(:,i1,ibond) 1614 !end do 1615 !end do 1616 !ENDDEBUG 1617 1618 iprim=1 1619 !first: bond values 1620 do ibond=1,deloc%nbond 1621 i1 = deloc%bonds(1,1,ibond) 1622 s1 = deloc%bonds(2,1,ibond) 1623 r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)& 1624 & +deloc%rshift(2,s1)*rprimd(:,2)& 1625 & +deloc%rshift(3,s1)*rprimd(:,3) 1626 i2 = deloc%bonds(1,2,ibond) 1627 s2 = deloc%bonds(2,2,ibond) 1628 r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)& 1629 & +deloc%rshift(2,s2)*rprimd(:,2)& 1630 & +deloc%rshift(3,s2)*rprimd(:,3) 1631 prim_int(iprim) = bond_length(r1,r2) 1632 iprim=iprim+1 1633 end do 1634 1635 !second: angle values (ang) 1636 do iang=1,deloc%nang 1637 i1 = deloc%angs(1,1,iang) 1638 s1 = deloc%angs(2,1,iang) 1639 r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)& 1640 & +deloc%rshift(2,s1)*rprimd(:,2)& 1641 & +deloc%rshift(3,s1)*rprimd(:,3) 1642 i2 = deloc%angs(1,2,iang) 1643 s2 = deloc%angs(2,2,iang) 1644 r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)& 1645 & +deloc%rshift(2,s2)*rprimd(:,2)& 1646 & +deloc%rshift(3,s2)*rprimd(:,3) 1647 i3 = deloc%angs(1,3,iang) 1648 s3 = deloc%angs(2,3,iang) 1649 r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)& 1650 & +deloc%rshift(2,s3)*rprimd(:,2)& 1651 & +deloc%rshift(3,s3)*rprimd(:,3) 1652 prim_int(iprim) = angle_ang(r1,r2,r3) 1653 iprim=iprim+1 1654 end do 1655 1656 !third: dihedral values 1657 do idihed=1,deloc%ndihed 1658 i1 = deloc%dihedrals(1,1,idihed) 1659 s1 = deloc%dihedrals(2,1,idihed) 1660 r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)& 1661 & +deloc%rshift(2,s1)*rprimd(:,2)& 1662 & +deloc%rshift(3,s1)*rprimd(:,3) 1663 i2 = deloc%dihedrals(1,2,idihed) 1664 s2 = deloc%dihedrals(2,2,idihed) 1665 r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)& 1666 & +deloc%rshift(2,s2)*rprimd(:,2)& 1667 & +deloc%rshift(3,s2)*rprimd(:,3) 1668 i3 = deloc%dihedrals(1,3,idihed) 1669 s3 = deloc%dihedrals(2,3,idihed) 1670 r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)& 1671 & +deloc%rshift(2,s3)*rprimd(:,2)& 1672 & +deloc%rshift(3,s3)*rprimd(:,3) 1673 i4 = deloc%dihedrals(1,4,idihed) 1674 s4 = deloc%dihedrals(2,4,idihed) 1675 r4(:) = xcart(:,i4)+deloc%rshift(1,s4)*rprimd(:,1)& 1676 & +deloc%rshift(2,s4)*rprimd(:,2)& 1677 & +deloc%rshift(3,s4)*rprimd(:,3) 1678 prim_int(iprim) = angle_dihedral(r1,r2,r3,r4) 1679 iprim=iprim+1 1680 end do 1681 1682 do icart=1,deloc%ncart 1683 prim_int(iprim) = xcart(deloc%carts(1,icart),deloc%carts(2,icart)) 1684 iprim=iprim+1 1685 end do 1686 1687 !DEBUG 1688 !write(std_out,*) 'Primitive internal coordinate values:' 1689 !do iprim=1,ninternal 1690 !if (iprim <= deloc%nbond) then 1691 !write(std_out,*) iprim, prim_int(iprim) 1692 !else 1693 !write(std_out,*) iprim, prim_int(iprim), prim_int(iprim)/pi*180.0_dp 1694 !end if 1695 !end do 1696 !ENDDEBUG 1697 1698 end subroutine calc_prim_int
m_abimover/delocint [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
delocint
FUNCTION
Datatype with the important variables in pred_delocint
NOTES
deloc <type(delocint)>=Important variables for pred_delocint | ! icenter = Index of the center of the number of shifts | nang = Number of angles | nbond = Number of bonds | ncart = Number of cartesian directions (used for constraints) | ndihed = Number of dihedrals | nrshift = Dimension of rshift | ninternal= Number of internal coordinates | ninternal=nbond+nang+ndihed+ncart | angs(2,3,nang) = Indexes to characterize angles | bonds(2,2,nbond)= For a bond between iatom and jatom | bonds(1,1,nbond) = iatom | bonds(2,1,nbond) = icenter | bonds(1,2,nbond) = jatom | bonds(2,2,nbond) = irshift | carts(2,ncart) = Index of total primitive internal, and atom (carts(2,:)) | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals | rshift(3,nrshift)= Shift in xred that must be done to find | all neighbors of a given atom within a | given number of neighboring shells
SOURCE
210 type,public :: delocint 211 212 ! scalars 213 integer :: icenter 214 integer :: nang 215 integer :: nbond 216 integer :: ncart 217 integer :: ndihed 218 integer :: nrshift 219 integer :: ninternal 220 221 ! arrays 222 integer,allocatable :: angs(:,:,:) 223 integer,allocatable :: bonds(:,:,:) 224 integer,allocatable :: carts(:,:) 225 integer,allocatable :: dihedrals(:,:,:) 226 real(dp),allocatable :: rshift(:,:) 227 228 end type delocint
m_abimover/delocint_fin [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
delocint_fin
FUNCTION
destructor function for delocint object INPUT deloc= container object for delocalized internal coordinates
OUTPUT
SOURCE
2274 subroutine delocint_fin(deloc) 2275 2276 type(delocint), intent(inout) :: deloc 2277 2278 ABI_SFREE(deloc%angs) 2279 ABI_SFREE(deloc%bonds) 2280 ABI_SFREE(deloc%carts) 2281 ABI_SFREE(deloc%dihedrals) 2282 ABI_SFREE(deloc%rshift) 2283 2284 end subroutine delocint_fin
m_abimover/delocint_ini [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
delocint_ini
FUNCTION
ini function for delocint object INPUT
OUTPUT
SIDE EFFECTS
deloc= container object for delocalized internal coordinates
SOURCE
2227 subroutine delocint_ini(deloc) 2228 2229 !Arguments ------------------------------------ 2230 !scalars 2231 type(delocint), intent(out) :: deloc 2232 2233 !Local variables ------------------------------ 2234 !scalars 2235 integer :: ii,irshift,jj,kk,nshell 2236 2237 ! ********************************************************************* 2238 2239 nshell=3 2240 deloc%nrshift=(2*nshell+1)**3 2241 deloc%icenter = nshell*(2*nshell+1)**2 + nshell*(2*nshell+1) + nshell + 1 2242 2243 ABI_MALLOC(deloc%rshift,(3,deloc%nrshift)) 2244 irshift=0 2245 do ii=-nshell,nshell 2246 do jj=-nshell,nshell 2247 do kk=-nshell,nshell 2248 irshift=irshift+1 2249 deloc%rshift(:,irshift) = (/dble(ii),dble(jj),dble(kk)/) 2250 end do 2251 end do 2252 end do 2253 2254 end subroutine delocint_ini
m_abimover/go_angles [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
go_angles
FUNCTION
Datatype all the information relevant to create angles between atoms inside and outside the cell
NOTES
This type is not used
SOURCE
379 type, public :: go_angles 380 381 !scalar 382 integer :: nangles ! Total number of bonds for the system 383 384 !arrays 385 integer,allocatable :: angle_vertex(:) ! Indices of the vertex atom 386 real(dp),allocatable :: angle_value(:) ! Value of angle in radians 387 real(dp),allocatable :: angle_bonds(:,:) ! Indices of the bonds 388 real(dp),allocatable :: angle_vect(:,:) ! Unitary vector perpendicular to the plane 389 390 end type go_angles 391 392 !public :: make_angles_new ! This routine is broken and should be tested before use.
m_abimover/go_bonds [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
go_bonds
FUNCTION
Datatype all the information relevant to create bonds between atoms inside and outside the cell
SOURCE
342 type, public :: go_bonds 343 344 !scalar 345 real(dp) :: tolerance ! To decide if consider bond the atom or not 346 ! 1.0 means that only consider values lower 347 ! than the sum of covalent radius 348 349 integer :: nbonds ! Total number of bonds for the system 350 351 !arrays 352 353 integer,allocatable :: nbondi(:) ! Number of bonds for atom i 354 integer,allocatable :: indexi(:,:) ! Indices of bonds for atom i 355 ! Positive: Vector from i to j 356 ! Negative: Vector from j to i 357 358 real(dp),allocatable :: bond_length(:) ! Bond lengths 359 real(dp),allocatable :: bond_vect(:,:) ! Unitary vectors for bonds 360 361 end type go_bonds
m_abimover/make_angles [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
make_angles
FUNCTION
INPUTS
OUTPUT
SOURCE
1240 subroutine make_angles(deloc,natom) 1241 1242 !Arguments ------------------------------------ 1243 !scalars 1244 integer,intent(in) :: natom 1245 type(delocint),intent(inout) :: deloc 1246 !arrays 1247 1248 !Local variables------------------------------- 1249 !scalars 1250 integer :: ia1,ia2,iang,ibond,is1,is2,ishift,ja1,ja2 1251 integer :: jbond,js1,js2 1252 !arrays 1253 integer,allocatable :: angs_tmp(:,:,:) 1254 1255 ! ************************************************************************* 1256 1257 !tentative first allocation: < 6 angles per bond. 1258 ABI_MALLOC(angs_tmp,(2,3,72*natom)) 1259 1260 deloc%nang = 0 1261 1262 do ibond=1, deloc%nbond 1263 ia1 = deloc%bonds(1,1,ibond) 1264 is1 = deloc%bonds(2,1,ibond) 1265 ia2 = deloc%bonds(1,2,ibond) 1266 is2 = deloc%bonds(2,2,ibond) 1267 do jbond=ibond+1,deloc%nbond 1268 ja1 = deloc%bonds(1,1,jbond) 1269 ja2 = deloc%bonds(1,2,jbond) 1270 do ishift=-(deloc%icenter-1),+(deloc%icenter-1) 1271 js1 = deloc%bonds(2,1,jbond)+ishift 1272 js2 = deloc%bonds(2,2,jbond)+ishift 1273 1274 if (ia1==ja1 .and. is1==js1) then 1275 deloc%nang = deloc%nang+1 1276 angs_tmp(:,1,deloc%nang) = (/ia2,is2/) 1277 angs_tmp(:,2,deloc%nang) = (/ia1,is1/) 1278 angs_tmp(:,3,deloc%nang) = (/ja2,js2/) 1279 1280 else if (ia1==ja2 .and. is1==js2) then 1281 deloc%nang = deloc%nang+1 1282 angs_tmp(:,1,deloc%nang) = (/ia2,is2/) 1283 angs_tmp(:,2,deloc%nang) = (/ia1,is1/) 1284 angs_tmp(:,3,deloc%nang) = (/ja1,js1/) 1285 1286 else if (ia2==ja2 .and. is2==js2) then 1287 deloc%nang = deloc%nang+1 1288 angs_tmp(:,1,deloc%nang) = (/ia1,is1/) 1289 angs_tmp(:,2,deloc%nang) = (/ia2,is2/) 1290 angs_tmp(:,3,deloc%nang) = (/ja1,js1/) 1291 1292 else if (ia2==ja1 .and. is2==js1) then 1293 deloc%nang = deloc%nang+1 1294 angs_tmp(:,1,deloc%nang) = (/ia1,is1/) 1295 angs_tmp(:,2,deloc%nang) = (/ia2,is2/) 1296 angs_tmp(:,3,deloc%nang) = (/ja2,js2/) 1297 1298 end if 1299 if (deloc%nang > 72*natom) then 1300 ABI_ERROR('too many angles found > 72*natom') 1301 end if 1302 end do 1303 end do ! jbond do 1304 end do ! ibond 1305 1306 ABI_SFREE(deloc%angs) 1307 ABI_MALLOC(deloc%angs,(2,3,deloc%nang)) 1308 do iang=1,deloc%nang 1309 deloc%angs(:,:,iang) = angs_tmp(:,:,iang) 1310 end do 1311 ABI_FREE(angs_tmp) 1312 1313 end subroutine make_angles
m_abimover/make_angles_new [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
make_angles_new
FUNCTION
Fill the contents of the angles structure, that contains all non redundant angles that could be generated between all the atoms in the unitary cell and their adjacent cells An angle is establish when an atom has two or more bonds. The angles structure contains information about the atoms involved, the value of the angle in radians, and the unitary vector perpendicular to the plane of the three atoms that build the angle.
INPUTS
natom= Number of atoms
ntypat= Number of type of atoms
rprimd= Dimensional primitive vectors of the cell
xcart= Cartesian coordinates of the atoms
znucl= Z number of the atom
bonds= Structure that store all the information about
bonds created by this routine:
nbonds= Total number of bonds
nbondi= Number of bonds for atom i
indexi= Indeces of bonds for atom i
bond_length= Distances between atoms i and j (including shift)
bond_vect= Unitary vector for the bond from i to j
tolerance= The tolerance is multiplied to the
adition of covalent radius to decide if a bond is created
OUTPUT
angles= Structure that store the information about
angles created by this routine
nangles= Total number of angles
angle_vertex= Index of the atom for that angle
angle_value= Value of the angle in radians
angle_bonds= Indices of the bonds
angle_vect= Unitary vector perpendicular to the plane of the angle
SOURCE
2330 !This routine has been disables since it's broken 2331 #if 0 2332 2333 subroutine make_angles_new(angles,bonds,natom,ntypat,rprimd,typat,xcart,znucl) 2334 2335 !Arguments ------------------------------------ 2336 !scalars 2337 integer,intent(in) :: natom,ntypat 2338 !arrays 2339 integer,intent(in) :: typat(natom) 2340 real(dp),intent(in) :: znucl(ntypat) 2341 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 2342 type(go_bonds),intent(in) :: bonds 2343 type(go_angles),intent(inout) :: angles 2344 2345 !Local variables ------------------------------ 2346 !scalars 2347 integer :: ii,jj,kk,iangle 2348 type(atomdata_t) :: atom 2349 2350 !arrays 2351 type(go_bonds) :: bonds_tmp 2352 character(len=2) :: symbol(ntypat) 2353 real(dp) :: amu(ntypat) 2354 integer :: shift(3,13) ! Represent all shift vectors that are not equivalent by central symmetry 2355 ! For example (1,1,-1) is equivalent to (-1,-1,1) 2356 ! It means that bond between atom i in the original cell and atom j in the 2357 ! cell with cordinates (1,1,-1) is equivalent to the bond between atom j in 2358 ! the orignal cell and atom i in the cell with coordinates (-1,-1,1) 2359 ! The trivial shift (0,0,0) is excluded here 2360 real(dp) :: rcov(ntypat) ! Covalent radius 2361 real(dp) :: rpt(3) 2362 2363 !*************************************************************************** 2364 !Beginning of executable session 2365 !*************************************************************************** 2366 ABI_ERROR("This routine is not tested") 2367 2368 !write(std_out,*) 'make_bonds 01' 2369 !########################################################## 2370 !### 01. Compute covalent radius 2371 2372 do ii=1,ntypat 2373 call atomdata_from_znucl(atom,znucl(ii)) 2374 amu(ii) = atom%amu 2375 rcov(ii) = atom%rcov 2376 symbol(ii) = symbol(ii) 2377 end do 2378 2379 !write(std_out,*) 'make_bonds 02' 2380 !########################################################## 2381 !### 02. Fill the 13 posible shift conecting adjacent cells 2382 2383 shift(:,:)=reshape( (/ 1,0,0,& 2384 & 0, 1, 0,& 2385 & 0, 0, 1,& 2386 & 1, 1, 0,& 2387 & 1,-1, 0,& 2388 & 0, 1, 1,& 2389 & 0, 1,-1,& 2390 & 1, 0, 1,& 2391 & 1, 0,-1,& 2392 & 1, 1, 1,& 2393 & 1,-1, 1,& 2394 & 1, 1,-1,& 2395 & 1,-1,-1 /), (/ 3, 13 /)) 2396 2397 !write(std_out,*) 'make_bonds 03' 2398 !########################################################## 2399 !### 03. Initialize the values of bonds 2400 2401 !The total number of bonds could not be predicted without 2402 !compute all the distances, but the extreme case is linking 2403 !all the atoms within all adjacent cells (natom*natom*13) 2404 !plus the all the bonds inside the original cell (natom*(natom-1)) 2405 2406 bonds_tmp%nbonds=0 2407 bonds_tmp%tolerance=bonds%tolerance 2408 ibond=0 2409 2410 ABI_MALLOC(bonds_tmp%bond_vect,(3,natom*natom*14-natom)) 2411 ABI_MALLOC(bonds_tmp%bond_length,(natom*natom*14-natom)) 2412 2413 !indexi contains the indices to the bonds 2414 ABI_MALLOC(bonds_tmp%indexi,(natom,natom*natom*14-natom)) 2415 2416 ABI_MALLOC(bonds_tmp%nbondi,(natom)) 2417 2418 bonds_tmp%indexi(:,:)=0 2419 bonds_tmp%nbondi(:)=0 2420 2421 !write(std_out,*) 'make_bonds 04' 2422 !########################################################## 2423 !### 04. Compute the bonds inside the original cell 2424 !### shift=(0,0,0) 2425 2426 do ii=1,natom 2427 rcov1 = rcov(typat(ii)) 2428 2429 do jj=ii+1,natom 2430 rcov2 = rcov(typat(jj)) 2431 2432 bl=bond_length(xcart(:,ii),xcart(:,jj)) 2433 2434 if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then 2435 ! We have a new bond, nbonds starts from 2436 ! 0, so it could be used to index the 2437 ! locations of bondij and distij 2438 2439 ! Increase the number of bonds 2440 bonds_tmp%nbonds= bonds_tmp%nbonds+1 2441 2442 ! The number of bonds for atoms ii and jj 2443 ! needs to raise by one 2444 bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1 2445 bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1 2446 2447 bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds 2448 ! The value for jj is negative to indicate that 2449 ! the vector is from ii to jj 2450 bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds 2451 2452 ! The unitary vector is always from ii to jj 2453 bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(xcart(:,jj)-xcart(:,ii))/bl 2454 bonds_tmp%bond_length(bonds_tmp%nbonds)=bl 2455 2456 end if 2457 2458 end do !! jj 2459 end do !! ii 2460 2461 !write(std_out,*) 'make_bonds 05' 2462 !########################################################## 2463 !### 05. Compute the bonds outside the original cell 2464 !### 13 shifts considered 2465 2466 !Bonds between identical atoms but in diferent cells are 2467 !allowed 2468 2469 do ii=1,natom 2470 rcov1 = rcov(typat(ii)) 2471 do jj=1,natom 2472 rcov2 = rcov(typat(jj)) 2473 2474 do irshift=1,13 2475 2476 do kk=1,3 2477 rpt(kk) = xcart(kk,jj)+& 2478 & shift(1,irshift)*rprimd(kk,1)+ & 2479 & shift(2,irshift)*rprimd(kk,2)+ & 2480 & shift(3,irshift)*rprimd(kk,3) 2481 end do 2482 2483 2484 bl =bond_length(xcart(:,ii),rpt) 2485 2486 if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then 2487 2488 ! We have a new bond, nbonds starts from 2489 ! 0, so it could be used to index the 2490 ! locations of bondij and distij 2491 2492 ! Increase the number of bonds 2493 bonds_tmp%nbonds= bonds_tmp%nbonds+1 2494 2495 ! The number of bonds for atoms ii and jj 2496 ! needs to raise by one 2497 bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1 2498 bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds 2499 2500 ! The value for jj is negative to indicate that 2501 ! the vector is from ii to jj 2502 bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1 2503 bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds 2504 2505 ! The unitary vector is always from ii to jj 2506 bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(rpt(:)-xcart(:,ii))/bl 2507 bonds_tmp%bond_length(bonds_tmp%nbonds)=bl 2508 2509 if (ii==jj) then 2510 bonds_tmp%nbonds= bonds_tmp%nbonds+1 2511 end if 2512 2513 end if 2514 2515 end do !! irshift 2516 2517 end do !! jj 2518 end do !! ii 2519 2520 call print_bonds(amu,bonds_tmp,natom,ntypat,symbol,typat,znucl) 2521 2522 2523 !write(std_out,*) 'make_bonds 05' 2524 !########################################################## 2525 !### 05. Deallocate all the arrays inside bonds 2526 !### allocate them with the right size and fill them 2527 2528 call bonds_free(bonds) 2529 2530 bonds%nbonds=bonds_tmp%nbonds 2531 2532 if (bonds%nbonds>0) then 2533 ! Allocate the arrays with exactly the rigth nbonds 2534 ABI_MALLOC(bonds%bond_vect,(3,bonds%nbonds)) 2535 ABI_MALLOC(bonds%bond_length,(bonds%nbonds)) 2536 ABI_MALLOC(bonds%indexi,(natom,bonds%nbonds)) 2537 ABI_MALLOC(bonds%nbondi,(natom)) 2538 2539 ! Fill the values 2540 bonds%bond_vect(:,1:bonds%nbonds)=bonds_tmp%bond_vect(:,1:bonds%nbonds) 2541 bonds%bond_length(1:bonds%nbonds)=bonds_tmp%bond_length(1:bonds%nbonds) 2542 bonds%indexi(:,1:bonds%nbonds)=bonds_tmp%indexi(:,1:bonds%nbonds) 2543 bonds%nbondi(:)=bonds_tmp%nbondi(:) 2544 end if 2545 2546 call bonds_free(bonds_tmp) 2547 2548 end subroutine make_angles_new
m_abimover/make_bonds [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
make_bonds
FUNCTION
INPUTS
OUTPUT
SOURCE
1464 subroutine make_bonds(deloc,natom,ntypat,rprimd,typat,xcart,znucl) 1465 1466 !Arguments ------------------------------------ 1467 !scalars 1468 integer,intent(in) :: natom,ntypat 1469 type(delocint),intent(inout) :: deloc 1470 !arrays 1471 integer,intent(in) :: typat(natom) 1472 real(dp),intent(in) :: znucl(:) ! znucl(ntypat) or 1473 ! znucl(npsp) ? 1474 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 1475 1476 !Local variables ------------------------------ 1477 !scalars 1478 integer :: iatom,ibond,irshift,itypat,jatom 1479 real(dp) :: bl,bondfudge,rcov1,rcov2 1480 type(atomdata_t) :: atom 1481 !arrays 1482 integer,allocatable :: bonds_tmp(:,:,:) 1483 real(dp) :: rcov(ntypat),rpt(3) 1484 1485 !************************************************************************ 1486 1487 do itypat=1,ntypat 1488 call atomdata_from_znucl(atom,znucl(itypat)) 1489 rcov(itypat) = atom%rcov 1490 end do 1491 1492 !write(std_out,*) ' rcov =', rcov 1493 !write(std_out,*) ' nrshift =', deloc%nrshift 1494 !write(std_out,*) ' xcart =', xcart 1495 !write(std_out,*) ' natom =',natom 1496 1497 !tentative first allocation: < 12 bonds per atom. 1498 ABI_MALLOC(bonds_tmp,(2,2,12*natom)) 1499 1500 bondfudge = 1.1_dp 1501 1502 deloc%nbond = 0 1503 1504 do iatom=1,natom 1505 rcov1 = rcov(typat(iatom)) 1506 do jatom=iatom+1,natom 1507 rcov2 = rcov(typat(jatom)) 1508 do irshift=1,deloc%nrshift 1509 rpt(:) = xcart(:,jatom) & 1510 & + deloc%rshift(1,irshift)*rprimd(:,1) & 1511 & + deloc%rshift(2,irshift)*rprimd(:,2) & 1512 & + deloc%rshift(3,irshift)*rprimd(:,3) 1513 bl = bond_length(xcart(:,iatom),rpt) 1514 1515 !write(std_out,*) ' bl, bondfudge*(rcov1+rcov2) = ',bl, bondfudge*(rcov1+rcov2) 1516 1517 if (bondfudge*(rcov1+rcov2) - bl > tol6) then 1518 deloc%nbond = deloc%nbond+1 1519 if (deloc%nbond > 12*natom) then 1520 ABI_ERROR('make_bonds: error too many bonds !') 1521 end if 1522 bonds_tmp(1,1,deloc%nbond) = iatom 1523 bonds_tmp(2,1,deloc%nbond) = deloc%icenter 1524 bonds_tmp(1,2,deloc%nbond) = jatom 1525 bonds_tmp(2,2,deloc%nbond) = irshift 1526 1527 !write(std_out,*) ' ibond bonds = ', deloc%nbond, bonds_tmp(:,:,deloc%nbond),xcart(:,iatom),rpt 1528 end if 1529 end do ! jatom 1530 end do 1531 end do ! iatom 1532 1533 ABI_SFREE(deloc%bonds) 1534 1535 ABI_MALLOC(deloc%bonds,(2,2,deloc%nbond)) 1536 do ibond=1,deloc%nbond 1537 deloc%bonds(:,:,ibond) = bonds_tmp(:,:,ibond) 1538 end do 1539 1540 ! do ibond=1,deloc%nbond 1541 ! write(std_out,*) ' make_bonds : bonds_tmp ', ibond, bonds_tmp(:,:,ibond) 1542 ! write(std_out,*) ' make_bonds : deloc%bonds ', ibond, deloc%bonds(:,:,ibond) 1543 ! end do 1544 1545 ABI_FREE(bonds_tmp) 1546 1547 end subroutine make_bonds
m_abimover/make_bonds_new [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
make_bonds_new
FUNCTION
Fill the contents of the bonds structure, that contains all non redundant bonds that could be generated between all the atoms in the unitary cell and their adjacent cells
INPUTS
natom= Number of atoms ntypat= Number of type of atoms rprimd= Dimensional primitive vectors of the cell xcart= Cartesian coordinates of the atoms znucl= Z number of the atom
OUTPUT
bonds= Structure that store all the information about
bonds created by this routine:
nbonds= Total number of bonds
nbondi= Number of bonds for atom i
indexi= Indeces of bonds for atom i
bond_length= Distances between atoms i and j (including shift)
bond_vect= Unitary vector for the bond from i to j
tolerance= The tolerance is multiplied to the
adition of covalent radius to decide if a bond is created
SOURCE
1886 subroutine make_bonds_new(bonds,natom,ntypat,rprimd,typat,xcart,znucl) 1887 1888 !Arguments ------------------------------------ 1889 !scalars 1890 integer,intent(in) :: natom,ntypat 1891 !arrays 1892 integer,intent(in) :: typat(natom) 1893 real(dp),intent(in) :: znucl(ntypat) 1894 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 1895 type(go_bonds),intent(inout) :: bonds 1896 1897 !Local variables ------------------------------ 1898 !scalars 1899 integer :: ii,jj,kk,ibond,irshift 1900 real(dp) :: rcov1,rcov2 1901 real(dp) :: bl 1902 type(go_bonds) :: bonds_tmp 1903 type(atomdata_t) :: atom 1904 1905 !arrays 1906 character(len=2) :: symbol(ntypat) 1907 real(dp) :: amu(ntypat) 1908 integer :: shift(3,13) ! Represent all shift vectors that are not equivalent by central symmetry 1909 ! For example (1,1,-1) is equivalent to (-1,-1,1) 1910 ! It means that bond between atom i in the original cell and atom j in the 1911 ! cell with cordinates (1,1,-1) is equivalent to the bond between atom j in 1912 ! the orignal cell and atom i in the cell with coordinates (-1,-1,1) 1913 ! The trivial shift (0,0,0) is excluded here 1914 real(dp) :: rcov(ntypat) ! Covalent radius 1915 real(dp) :: rpt(3) 1916 1917 !*************************************************************************** 1918 !Beginning of executable session 1919 !*************************************************************************** 1920 1921 !write(std_out,*) 'make_bonds 01' 1922 !########################################################## 1923 !### 01. Compute covalent radius 1924 1925 do ii=1,ntypat 1926 call atomdata_from_znucl(atom,znucl(ii)) 1927 amu(ii) = atom%amu 1928 rcov(ii) = atom%rcov 1929 symbol(ii) = atom%symbol 1930 end do 1931 1932 !write(std_out,*) 'make_bonds 02' 1933 !########################################################## 1934 !### 02. Fill the 13 posible shift conecting adjacent cells 1935 1936 shift(:,:)=reshape( (/ 1,0,0,& 1937 & 0, 1, 0,& 1938 & 0, 0, 1,& 1939 & 1, 1, 0,& 1940 & 1,-1, 0,& 1941 & 0, 1, 1,& 1942 & 0, 1,-1,& 1943 & 1, 0, 1,& 1944 & 1, 0,-1,& 1945 & 1, 1, 1,& 1946 & 1,-1, 1,& 1947 & 1, 1,-1,& 1948 & 1,-1,-1 /), (/ 3, 13 /)) 1949 1950 !write(std_out,*) 'make_bonds 03' 1951 !########################################################## 1952 !### 03. Initialize the values of bonds 1953 1954 !The total number of bonds could not be predicted without 1955 !compute all the distances, but the extreme case is linking 1956 !all the atoms within all adjacent cells (natom*natom*13) 1957 !plus the all the bonds inside the original cell (natom*(natom-1)) 1958 1959 bonds_tmp%nbonds=0 1960 bonds_tmp%tolerance=bonds%tolerance 1961 ibond=0 1962 1963 ABI_MALLOC(bonds_tmp%bond_vect,(3,natom*natom*14-natom)) 1964 ABI_MALLOC(bonds_tmp%bond_length,(natom*natom*14-natom)) 1965 1966 !indexi contains the indices to the bonds 1967 ABI_MALLOC(bonds_tmp%indexi,(natom,natom*natom*14-natom)) 1968 1969 ABI_MALLOC(bonds_tmp%nbondi,(natom)) 1970 1971 bonds_tmp%indexi(:,:)=0 1972 bonds_tmp%nbondi(:)=0 1973 bonds_tmp%bond_vect(:,:)=0.0 1974 bonds_tmp%bond_length(:)=0.0 1975 1976 !write(std_out,*) 'make_bonds 04' 1977 !########################################################## 1978 !### 04. Compute the bonds inside the original cell 1979 !### shift=(0,0,0) 1980 1981 do ii=1,natom 1982 rcov1 = rcov(typat(ii)) 1983 1984 do jj=ii+1,natom 1985 rcov2 = rcov(typat(jj)) 1986 1987 bl=bond_length(xcart(:,ii),xcart(:,jj)) 1988 1989 if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then 1990 ! We have a new bond, nbonds starts from 1991 ! 0, so it could be used to index the 1992 ! locations of bondij and distij 1993 1994 ! Increase the number of bonds 1995 bonds_tmp%nbonds= bonds_tmp%nbonds+1 1996 1997 ! The number of bonds for atoms ii and jj 1998 ! needs to raise by one 1999 bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1 2000 bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1 2001 2002 bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds 2003 ! The value for jj is negative to indicate that 2004 ! the vector is from ii to jj 2005 bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds 2006 2007 ! The unitary vector is always from ii to jj 2008 bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(xcart(:,jj)-xcart(:,ii))/bl 2009 bonds_tmp%bond_length(bonds_tmp%nbonds)=bl 2010 2011 end if 2012 2013 end do !! jj 2014 end do !! ii 2015 2016 !write(std_out,*) 'make_bonds 05' 2017 !########################################################## 2018 !### 05. Compute the bonds outside the original cell 2019 !### 13 shifts considered 2020 2021 !Bonds between identical atoms but in diferent cells are 2022 !allowed 2023 2024 do ii=1,natom 2025 rcov1 = rcov(typat(ii)) 2026 do jj=1,natom 2027 rcov2 = rcov(typat(jj)) 2028 2029 do irshift=1,13 2030 2031 do kk=1,3 2032 rpt(kk) = xcart(kk,jj)+& 2033 & shift(1,irshift)*rprimd(kk,1)+ & 2034 & shift(2,irshift)*rprimd(kk,2)+ & 2035 & shift(3,irshift)*rprimd(kk,3) 2036 end do 2037 2038 2039 bl =bond_length(xcart(:,ii),rpt) 2040 2041 if (bonds_tmp%tolerance*(rcov1+rcov2) > bl) then 2042 2043 ! We have a new bond, nbonds starts from 2044 ! 0, so it could be used to index the 2045 ! locations of bondij and distij 2046 2047 ! Increase the number of bonds 2048 bonds_tmp%nbonds= bonds_tmp%nbonds+1 2049 2050 ! The number of bonds for atoms ii and jj 2051 ! needs to raise by one 2052 bonds_tmp%nbondi(ii)= bonds_tmp%nbondi(ii)+1 2053 bonds_tmp%indexi(ii,bonds_tmp%nbondi(ii))=bonds_tmp%nbonds 2054 2055 ! The value for jj is negative to indicate that 2056 ! the vector is from ii to jj 2057 bonds_tmp%nbondi(jj)= bonds_tmp%nbondi(jj)+1 2058 bonds_tmp%indexi(jj,bonds_tmp%nbondi(jj))=-bonds_tmp%nbonds 2059 2060 ! The unitary vector is always from ii to jj 2061 bonds_tmp%bond_vect(:,bonds_tmp%nbonds)=(rpt(:)-xcart(:,ii))/bl 2062 bonds_tmp%bond_length(bonds_tmp%nbonds)=bl 2063 2064 if (ii==jj) then 2065 bonds_tmp%nbonds= bonds_tmp%nbonds+1 2066 end if 2067 2068 end if 2069 2070 end do !! irshift 2071 2072 end do !! jj 2073 end do !! ii 2074 2075 call print_bonds(amu,bonds_tmp,natom,ntypat,symbol,typat,znucl) 2076 2077 2078 !write(std_out,*) 'make_bonds 05' 2079 !########################################################## 2080 !### 05. Deallocate all the arrays inside bonds 2081 !### allocate them with the right size and fill them 2082 2083 call bonds_free(bonds) 2084 2085 bonds%nbonds=bonds_tmp%nbonds 2086 2087 if (bonds%nbonds>0) then 2088 ! Allocate the arrays with exactly the rigth nbonds 2089 ABI_MALLOC(bonds%bond_vect,(3,bonds%nbonds)) 2090 ABI_MALLOC(bonds%bond_length,(bonds%nbonds)) 2091 ABI_MALLOC(bonds%indexi,(natom,bonds%nbonds)) 2092 ABI_MALLOC(bonds%nbondi,(natom)) 2093 2094 ! Fill the values 2095 bonds%bond_vect(:,1:bonds%nbonds)=bonds_tmp%bond_vect(:,1:bonds%nbonds) 2096 bonds%bond_length(1:bonds%nbonds)=bonds_tmp%bond_length(1:bonds%nbonds) 2097 bonds%indexi(:,1:bonds%nbonds)=bonds_tmp%indexi(:,1:bonds%nbonds) 2098 bonds%nbondi(:)=bonds_tmp%nbondi(:) 2099 end if 2100 2101 call bonds_free(bonds_tmp) 2102 2103 end subroutine make_bonds_new
m_abimover/make_dihedrals [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
make_dihedrals
FUNCTION
INPUTS
OUTPUT
SOURCE
1330 subroutine make_dihedrals(badangles,deloc) 1331 1332 !Arguments ------------------------------------ 1333 !scalars 1334 type(delocint),intent(inout) :: deloc 1335 !arrays 1336 integer,intent(in) :: badangles(deloc%nang) 1337 1338 !Local variables------------------------------- 1339 !scalars 1340 integer :: chkdihed,ia1,ia2,ia3,iang,idihed,is1,is2 1341 integer :: is3,ishift,ja1,ja2,ja3,jang,js1,js2,js3,maxshift 1342 integer :: minshift 1343 !arrays 1344 integer,allocatable :: diheds_tmp(:,:,:) 1345 1346 ! ************************************************************************* 1347 1348 !tentative first allocation: < 6 dihedrals per angle. 1349 ABI_MALLOC(diheds_tmp,(2,4,6*deloc%nang)) 1350 1351 deloc%ndihed = 0 1352 diheds_tmp(:,:,:) = 0 1353 1354 do iang=1,deloc%nang 1355 if (badangles(iang) == 1) cycle 1356 ia1 = deloc%angs(1,1,iang) 1357 is1 = deloc%angs(2,1,iang) 1358 ia2 = deloc%angs(1,2,iang) 1359 is2 = deloc%angs(2,2,iang) 1360 ia3 = deloc%angs(1,3,iang) 1361 is3 = deloc%angs(2,3,iang) 1362 1363 do jang=iang+1,deloc%nang 1364 if (badangles(jang) == 1) cycle 1365 ja1 = deloc%angs(1,1,jang) 1366 ja2 = deloc%angs(1,2,jang) 1367 ja3 = deloc%angs(1,3,jang) 1368 do ishift=-(deloc%icenter-1),(deloc%icenter-1) 1369 js1 = deloc%angs(2,1,jang)+ishift 1370 js2 = deloc%angs(2,2,jang)+ishift 1371 js3 = deloc%angs(2,3,jang)+ishift 1372 1373 chkdihed=0 1374 if (ia2==ja1 .and. is2==js1) then 1375 if (ia1==ja2 .and. is1==js2) then 1376 deloc%ndihed = deloc%ndihed+1 1377 diheds_tmp(:,1,deloc%ndihed) = (/ia3,is3/) 1378 diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/) 1379 diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/) 1380 diheds_tmp(:,4,deloc%ndihed) = (/ja3,js3/) 1381 chkdihed=1 1382 else if (ia3==ja2 .and. is3==js2) then 1383 deloc%ndihed = deloc%ndihed+1 1384 diheds_tmp(:,1,deloc%ndihed) = (/ia1,is1/) 1385 diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/) 1386 diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/) 1387 diheds_tmp(:,4,deloc%ndihed) = (/ja3,js3/) 1388 chkdihed=1 1389 end if 1390 else if (ia2==ja3 .and. is2==js3) then 1391 if (ia1==ja2 .and. is1==js2) then 1392 deloc%ndihed = deloc%ndihed+1 1393 diheds_tmp(:,1,deloc%ndihed) = (/ia3,is3/) 1394 diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/) 1395 diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/) 1396 diheds_tmp(:,4,deloc%ndihed) = (/ja1,js1/) 1397 chkdihed=1 1398 else if (ia3==ja2 .and. is3==js2) then 1399 deloc%ndihed = deloc%ndihed+1 1400 diheds_tmp(:,1,deloc%ndihed) = (/ia1,is1/) 1401 diheds_tmp(:,2,deloc%ndihed) = (/ia2,is2/) 1402 diheds_tmp(:,3,deloc%ndihed) = (/ja2,js2/) 1403 diheds_tmp(:,4,deloc%ndihed) = (/ja1,js1/) 1404 chkdihed=1 1405 end if 1406 end if 1407 if (deloc%ndihed > 6*deloc%nang) then 1408 ABI_ERROR('make_dihedrals : too many dihedrals found > 6*nang') 1409 end if 1410 if (chkdihed == 1) then 1411 if ( diheds_tmp(1,4,deloc%ndihed) == diheds_tmp(1,1,deloc%ndihed) .and.& 1412 & diheds_tmp(2,4,deloc%ndihed) == diheds_tmp(2,1,deloc%ndihed) ) then 1413 write(std_out,*) 'make_dihedrals : Bad dihedral was found: atom1 == atom4. Discarding.' 1414 diheds_tmp(:,:,deloc%ndihed) = 0 1415 deloc%ndihed = deloc%ndihed-1 1416 end if 1417 end if 1418 end do 1419 end do 1420 ! end jang do 1421 end do 1422 !end iang do 1423 1424 ABI_SFREE(deloc%dihedrals) 1425 1426 ABI_MALLOC(deloc%dihedrals,(2,4,deloc%ndihed)) 1427 do idihed=1,deloc%ndihed 1428 deloc%dihedrals(:,:,idihed) = diheds_tmp(:,:,idihed) 1429 1430 ! minshift = minval(diheds_tmp(2,:,idihed)) 1431 ! if (minshift <= 0) then 1432 ! deloc%dihedrals(2,:,idihed) = deloc%dihedrals(2,:,idihed)+minshift+1 1433 ! end if 1434 ! maxshift = maxval(diheds_tmp(2,:,idihed)) 1435 ! if (maxshift > deloc%nrshift) then 1436 ! deloc%dihedrals(2,:,idihed) = deloc%dihedrals(2,:,idihed)-maxshift 1437 ! end if 1438 ! 1439 minshift = minval(diheds_tmp(2,:,idihed)) 1440 maxshift = maxval(diheds_tmp(2,:,idihed)) 1441 if (minshift <= 0 .or. maxshift > deloc%nrshift) then 1442 ABI_ERROR("dihedral extends beyond first neighboring unit cells!") 1443 end if 1444 end do 1445 ABI_FREE(diheds_tmp) 1446 1447 end subroutine make_dihedrals
m_abimover/make_prim_internals [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
make_prim_internals
FUNCTION
Determine the bonds, angles and dihedrals for a starting geometry, based on covalent radii for the atoms.
INPUTS
natom = Number of atoms (dtset%natom) nrshift= dimension of rshift rprimd(3,3)=dimensional real space primitive translations (bohr) rshift(3,nrshift)=shift in xred that must be done to find all neighbors of a given atom within a given number of neighboring shells xcart(3,natom)=cartesian coordinates of atoms (bohr)
OUTPUT
SIDE EFFECTS
deloc <type(delocint)>=Important variables for | pred_delocint ! icenter = Index of the center of the number of shifts | nang = Number of angles | nbond = Number of bonds | ncart = Number of cartesian directions | (used for constraints) | ndihed = Number of dihedrals | nrshift = Dimension of rshift | ninternal= Number of internal coordinates | ninternal=nbond+nang+ndihed+ncart | | angs(2,3,nang) = Indexes to characterize angles | bonds(2,2,nbond)= For a bond between iatom and jatom | bonds(1,1,nbond) = iatom | bonds(2,1,nbond) = icenter | bonds(1,2,nbond) = jatom | bonds(2,2,nbond) = irshift | carts(2,ncart) = Index of total primitive internal, | and atom (carts(2,:)) | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals | | rshift(3,nrshift)= Shift in xred that must be done to find | all neighbors of a given atom within a | given number of neighboring shells
NOTES
Adds cartesian coordinates if the number of internals with a given atom is < 4 the chosen coordinate could be optimized to be less dependent of the internals already incorporated.
SOURCE
1085 subroutine make_prim_internals(deloc,natom,ntypat,rprimd,typat,xcart,znucl) 1086 1087 !Arguments ------------------------------------ 1088 !scalars 1089 type(delocint),intent(inout) :: deloc 1090 integer,intent(in) :: natom,ntypat 1091 !arrays 1092 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 1093 integer,intent(in) :: typat(natom) 1094 real(dp),intent(in) :: znucl(:) ! znucl(ntypat) or znucl(npsp) ? 1095 1096 !Local variables ------------------------------ 1097 ! function 1098 !scalars 1099 integer :: iang,iatom,ibond,icart,idihed,ii 1100 real(dp) :: spp 1101 !arrays 1102 integer :: particip_atom(natom) 1103 integer,allocatable :: badangles(:) 1104 real(dp) :: rpt1(3),rpt3(3) !rpt2(3) 1105 1106 !************************************************************************ 1107 1108 particip_atom(:) = 0 1109 1110 call make_bonds(deloc,natom,ntypat,rprimd,typat,xcart,znucl) 1111 1112 do ibond=1,deloc%nbond 1113 write(std_out,'(a,i4,2(2i5,2x))') 'bond ', ibond, deloc%bonds(:,:,ibond) 1114 particip_atom(deloc%bonds(1,1,ibond)) = particip_atom(deloc%bonds(1,1,ibond))+1 1115 particip_atom(deloc%bonds(1,2,ibond)) = particip_atom(deloc%bonds(1,2,ibond))+1 1116 end do 1117 1118 call make_angles(deloc,natom) 1119 1120 ABI_MALLOC(badangles,(deloc%nang)) 1121 badangles(:) = 0 1122 do iang=1,deloc%nang 1123 write(std_out,'(a,i4,3(2i5,2x))') 'angle ', iang, deloc%angs(:,:,iang) 1124 particip_atom(deloc%angs(1,1,iang)) = particip_atom(deloc%angs(1,1,iang))+1 1125 particip_atom(deloc%angs(1,2,iang)) = particip_atom(deloc%angs(1,2,iang))+1 1126 particip_atom(deloc%angs(1,3,iang)) = particip_atom(deloc%angs(1,3,iang))+1 1127 1128 ! DEBUG 1129 ! rpt1(:) = xcart(:,deloc%angs(1,1,iang)) & 1130 ! & + deloc%rshift(1,deloc%angs(2,1,iang))*rprimd(:,1) & 1131 ! & + deloc%rshift(2,deloc%angs(2,1,iang))*rprimd(:,2) & 1132 ! & + deloc%rshift(3,deloc%angs(2,1,iang))*rprimd(:,3) 1133 ! rpt2(:) = xcart(:,deloc%angs(1,2,iang)) & 1134 ! & + deloc%rshift(1,deloc%angs(2,2,iang))*rprimd(:,1) & 1135 ! & + deloc%rshift(2,deloc%angs(2,2,iang))*rprimd(:,2) & 1136 ! & + deloc%rshift(3,deloc%angs(2,2,iang))*rprimd(:,3) 1137 ! rpt3(:) = xcart(:,deloc%angs(1,3,iang)) & 1138 ! & + deloc%rshift(1,deloc%angs(2,3,iang))*rprimd(:,1) & 1139 ! & + deloc%rshift(2,deloc%angs(2,3,iang))*rprimd(:,2) & 1140 ! & + deloc%rshift(3,deloc%angs(2,3,iang))*rprimd(:,3) 1141 ! write(std_out,*) rpt1,rpt2,rpt3,bond_length(rpt1,rpt2),bond_length(rpt2,rpt3) 1142 ! ENDDEBUG 1143 1144 ! check if angles are 180 degrees: discard the dihedrals in that case. 1145 rpt1(:) = xcart(:,deloc%angs(1,1,iang)) & 1146 & + deloc%rshift(1,deloc%angs(2,1,iang))*rprimd(:,1) & 1147 & + deloc%rshift(2,deloc%angs(2,1,iang))*rprimd(:,2) & 1148 & + deloc%rshift(3,deloc%angs(2,1,iang))*rprimd(:,3) & 1149 & - xcart(:,deloc%angs(1,2,iang)) & 1150 & - deloc%rshift(1,deloc%angs(2,2,iang))*rprimd(:,1) & 1151 & - deloc%rshift(2,deloc%angs(2,2,iang))*rprimd(:,2) & 1152 & - deloc%rshift(3,deloc%angs(2,2,iang))*rprimd(:,3) 1153 1154 rpt3(:) = xcart(:,deloc%angs(1,3,iang)) & 1155 & + deloc%rshift(1,deloc%angs(2,3,iang))*rprimd(:,1) & 1156 & + deloc%rshift(2,deloc%angs(2,3,iang))*rprimd(:,2) & 1157 & + deloc%rshift(3,deloc%angs(2,3,iang))*rprimd(:,3) & 1158 & - xcart(:,deloc%angs(1,2,iang)) & 1159 & - deloc%rshift(1,deloc%angs(2,2,iang))*rprimd(:,1) & 1160 & - deloc%rshift(2,deloc%angs(2,2,iang))*rprimd(:,2) & 1161 & - deloc%rshift(3,deloc%angs(2,2,iang))*rprimd(:,3) 1162 spp = (rpt1(1)*rpt3(1)+rpt1(2)*rpt3(2)+rpt1(3)*rpt3(3))& 1163 & / sqrt(rpt1(1)*rpt1(1)+rpt1(2)*rpt1(2)+rpt1(3)*rpt1(3)) & 1164 & / sqrt(rpt3(1)*rpt3(1)+rpt3(2)*rpt3(2)+rpt3(3)*rpt3(3)) 1165 if (abs(abs(spp) - one) < tol6) then 1166 write(std_out,*) 'make_prim_internals : an angle is too close to 180 degrees:' 1167 write(std_out,*) ' will discard dihedrals using it ' 1168 badangles(iang) = 1 1169 end if 1170 end do 1171 1172 call make_dihedrals(badangles,deloc) 1173 ABI_FREE(badangles) 1174 1175 do idihed=1,deloc%ndihed 1176 write(std_out,'(a,i4,4(2i5,2x))') 'dihedral ', idihed, deloc%dihedrals(:,:,idihed) 1177 particip_atom(deloc%dihedrals(1,1,idihed)) = particip_atom(deloc%dihedrals(1,1,idihed))+1 1178 particip_atom(deloc%dihedrals(1,2,idihed)) = particip_atom(deloc%dihedrals(1,2,idihed))+1 1179 particip_atom(deloc%dihedrals(1,3,idihed)) = particip_atom(deloc%dihedrals(1,3,idihed))+1 1180 particip_atom(deloc%dihedrals(1,4,idihed)) = particip_atom(deloc%dihedrals(1,4,idihed))+1 1181 1182 ! do ii=1,4 1183 ! write(std_out,'((3E16.6,2x))') xcart(:,deloc%dihedrals(1,ii,idihed)) + & 1184 ! & deloc%rshift(1,deloc%dihedrals(2,ii,idihed))*rprimd(:,1) + & 1185 ! & deloc%rshift(2,deloc%dihedrals(2,ii,idihed))*rprimd(:,2) + & 1186 ! & deloc%rshift(2,deloc%dihedrals(2,ii,idihed))*rprimd(:,3) 1187 ! end do 1188 end do 1189 1190 write(std_out,*) 'make_deloc_internals: nbond,nang,ndihed = ', deloc%nbond,deloc%nang,deloc%ndihed 1191 1192 !Check all atoms participate in at least 4 primitives. Otherwise, we should 1193 !probably add cartesian coordinates to the internal ones. 1194 deloc%ncart = 0 1195 do iatom=1,natom 1196 if (particip_atom(iatom) < 4) then 1197 write(std_out,*) ' make_prim_internals : Warning : atom ', iatom, & 1198 & ' does not belong to enough primitives to determine its' 1199 write(std_out,*) ' position uniquely ! instead : ', particip_atom(iatom) 1200 write(std_out,*) ' Will add cartesian coordinates to set of internals.' 1201 ! write(std_out,*) ' Not done yet.' 1202 ! stop 1203 deloc%ncart = deloc%ncart + 4-particip_atom(iatom) 1204 end if 1205 end do 1206 ABI_SFREE(deloc%carts) 1207 ABI_MALLOC(deloc%carts ,(2,deloc%ncart)) 1208 icart = 0 1209 do iatom=1,natom 1210 if (particip_atom(iatom) < 4) then 1211 ! kind of arbitrary : include first few directions for the atom: x, then y then z 1212 do ii=1,4-particip_atom(iatom) 1213 icart = icart+1 1214 deloc%carts(1,icart) = ii 1215 deloc%carts(2,icart) = iatom 1216 end do 1217 end if 1218 end do 1219 1220 !ninternal=nbond+nang+ndihed 1221 deloc%ninternal=deloc%nbond+deloc%nang+deloc%ndihed+deloc%ncart 1222 1223 end subroutine make_prim_internals
m_abimover/mttk_fin [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
mttk_fin
FUNCTION
destructor function for mttk object INPUT mttk
OUTPUT
SOURCE
964 subroutine mttk_fin(mttk_vars) 965 966 type(mttk_type), intent(inout) :: mttk_vars 967 968 ABI_SFREE(mttk_vars%glogs) 969 ABI_SFREE(mttk_vars%vlogs) 970 ABI_SFREE(mttk_vars%xlogs) 971 972 end subroutine mttk_fin
m_abimover/mttk_ini [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
mttk_ini
FUNCTION
destructor function for mttk object INPUT mttk
OUTPUT
SOURCE
935 subroutine mttk_ini(mttk_vars,nnos) 936 937 integer,intent(in) :: nnos 938 type(mttk_type), intent(out) :: mttk_vars 939 940 ABI_MALLOC(mttk_vars%glogs,(nnos)) 941 ABI_MALLOC(mttk_vars%vlogs,(nnos)) 942 ABI_MALLOC(mttk_vars%xlogs,(nnos)) 943 944 end subroutine mttk_ini
m_abimover/mttk_type [ Types ]
[ Top ] [ m_abimover ] [ Types ]
NAME
mttk_type
FUNCTION
For Martyna et al. (TTK) reversible MD integration scheme and related data
SOURCE
242 type, public :: mttk_type 243 244 real(dp) :: glogv 245 !Logarithm of the volume 246 247 real(dp) :: vlogv 248 !Derivative of logv 249 250 real(dp) :: gboxg(3,3) 251 !Imbalance in pressure (see paper) 252 253 real(dp) :: vboxg(3,3) 254 !Velocity of log rprimd (see paper) 255 256 real(dp), allocatable :: glogs(:) 257 ! glogs(nnos) 258 ! Imbalance of kinetic energy 259 260 real(dp), allocatable :: vlogs(:) 261 ! vlogs(nnos) 262 ! Velocities of thermostat variables 263 264 real(dp), allocatable :: xlogs(:) 265 ! xlogs(nnos) 266 ! Positions of thermostat variables 267 268 end type mttk_type
m_abimover/print_bonds [ Functions ]
[ Top ] [ m_abimover ] [ Functions ]
NAME
print_bonds
FUNCTION
Print the bonds
INPUTS
natom= Number of atoms ntypat= Number of type of atoms znucl= Z number of the atom
OUTPUT
bonds= Structure that store all the information about
bonds created by this routine:
nbonds= Total number of bonds
bondij= Unitary vector along the bond direction
distij= Distances between atoms i and j (including shift)
listij= Indices of bonds going from i to j
listji= Indices of bonds going from j to i
indexij= Number of bonds between i and j
indexji= Number of bonds between j and i
tolerance
SOURCE
2160 subroutine print_bonds(amu,bonds,natom,ntypat,symbol,typat,znucl) 2161 2162 !Arguments ------------------------------------ 2163 !scalars 2164 integer,intent(in) :: natom,ntypat 2165 integer,intent(in) :: typat(natom) 2166 real(dp),intent(in) :: znucl(ntypat) 2167 real(dp),intent(in) :: amu(ntypat) 2168 character(len=2),intent(in) :: symbol(ntypat) 2169 type(go_bonds),intent(in) :: bonds 2170 2171 !Local variables ------------------------------ 2172 !scalars 2173 integer :: ii,jj,kk 2174 2175 ! ********************************************************************* 2176 2177 write(std_out,'(a)') ch10 2178 write(std_out,'(a,72a,a)') '---BONDS',('-',kk=1,72),ch10 2179 write(std_out,'(a,i3)') 'Number of atoms: ',natom 2180 write(std_out,'(a,i3)') 'Number of bonds: ',bonds%nbonds 2181 write(std_out,'(a,f6.3,a,a)') 'Tolerance of bonds: ',bonds%tolerance,' times the sum of covalent radius',ch10 2182 2183 do ii=1,natom 2184 write(std_out,'(a,i3)') 'ATOM number: ',ii 2185 write(std_out,'(a,f8.3)') ' Z: ',znucl(typat(ii)) 2186 write(std_out,'(a,f8.3)') ' Weight: ',amu(typat(ii)) 2187 write(std_out,'(a,a3)') ' Symbol: ',symbol(typat(ii)) 2188 write(std_out,'(a,i3)') ' Number of bonds: ',bonds%nbondi(ii) 2189 2190 do jj=1,bonds%nbondi(ii) 2191 write(std_out,'(a,i3,a,a,i3,a,3f7.3,a,f7.3)') ' [',jj,']',& 2192 & ' Index of bond: ',bonds%indexi(ii,jj),& 2193 & ' Unitary vector: ',bonds%bond_vect(:,abs(bonds%indexi(ii,jj))),& 2194 & ' Bond length: ',bonds%bond_length(abs(bonds%indexi(ii,jj))) 2195 end do 2196 2197 end do 2198 2199 do ii=1,bonds%nbonds 2200 write(std_out,'(a,i3)') 'BOND Index=',ii 2201 write(std_out,'(a,3f8.3)') ' Vector',bonds%bond_vect(:,ii) 2202 write(std_out,'(a,f8.3)') ' bond Length',bonds%bond_length(ii) 2203 end do 2204 2205 end subroutine print_bonds