diff --git a/Makefile b/Makefile index 3c675b1d..560509aa 100644 --- a/Makefile +++ b/Makefile @@ -33,6 +33,7 @@ #! HND X #! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +CUSTOM_F95FLAGS += -frealloc-lhs ifeq (${QUIP_ARCH},) include Makefile.arch @@ -54,7 +55,7 @@ SOAP_TURBO_F95_FILES = soap_turbo_functions soap_turbo_radial soap_turbo_angular SOAP_TURBO_F95_SOURCES = ${addsuffix .f90, ${SOAP_TURBO_F95_FILES}} SOAP_TURBO_F95_OBJS = ${addsuffix .o, ${SOAP_TURBO_F95_FILES}} -GAP1_F95_FILES += descriptors gp_predict descriptors_wrapper clustering +GAP1_F95_FILES += find_water_triplets_noncommercial descriptors gp_predict descriptors_wrapper clustering GAP1_F95_SOURCES = ${addsuffix .f95, ${GAP1_F95_FILES}} GAP1_F95_OBJS = ${addsuffix .o, ${GAP1_F95_FILES}} diff --git a/descriptors.f95 b/descriptors.f95 index f9d1028b..7917c65e 100644 --- a/descriptors.f95 +++ b/descriptors.f95 @@ -107,6 +107,7 @@ module descriptors_module integer, parameter, public :: DT_DISTANCE_NB = 29 integer, parameter, public :: DT_SOAP_EXPRESS = 30 integer, parameter, public :: DT_SOAP_TURBO = 31 + integer, parameter, public :: DT_WATER_TRIMER = 32 integer, parameter :: NP_WATER_DIMER = 8 integer, parameter :: NP_A2_DIMER = 8 @@ -485,7 +486,7 @@ module descriptors_module public :: soap, general_monomer, bispectrum_so4, bispectrum_so3, behler, distance_2b, & coordination, angle_3b, co_angle_3b, co_distance_2b, cosnx, trihis, water_monomer, & water_dimer, a2_dimer, bond_real_space, power_so3, power_so4, an_monomer, general_dimer, & - general_trimer, rdf, as_distance_2b, molecule_lo_d, alex, com_dimer, distance_nb, & + general_trimer, water_trimer, rdf, as_distance_2b, molecule_lo_d, alex, com_dimer, distance_nb, & descriptor_data_mono, fourier_so4_type, radialfunction_type, transfer_parameters_type, & ab_dimer, atom_real_space, spherical_harmonics_type, behler_g2, behler_g3, soap_turbo, soap_express #else @@ -533,6 +534,7 @@ module descriptors_module type(general_monomer) :: descriptor_general_monomer type(general_dimer) :: descriptor_general_dimer type(general_trimer) :: descriptor_general_trimer + type(water_trimer) :: descriptor_water_trimer type(molecule_lo_d) :: descriptor_molecule_lo_d type(com_dimer) :: descriptor_com_dimer type(soap_express) :: descriptor_soap_express @@ -577,7 +579,7 @@ module descriptors_module coordination_initialise, angle_3b_initialise, co_angle_3b_initialise, co_distance_2b_initialise, cosnx_initialise, trihis_initialise, & water_monomer_initialise, water_dimer_initialise, A2_dimer_initialise, AB_dimer_initialise, distance_Nb_initialise, rdf_initialise, as_distance_2b_initialise, alex_initialise, & atom_real_space_initialise, power_so3_initialise, power_SO4_initialise, soap_initialise, soap_turbo_initialise, & - general_monomer_initialise, general_dimer_initialise, general_trimer_initialise, molecule_lo_d_initialise, AN_monomer_initialise, & + general_monomer_initialise, general_dimer_initialise, general_trimer_initialise, water_trimer_initialise, molecule_lo_d_initialise, AN_monomer_initialise, & bond_real_space_initialise, transfer_initialise, com_dimer_initialise, soap_express_initialise #else module procedure descriptor_initialise, RadialFunction_initialise, fourier_so4_initialise, & @@ -595,7 +597,7 @@ module descriptors_module bispectrum_SO4_finalise, bispectrum_SO3_finalise, behler_finalise, distance_2b_finalise, coordination_finalise, angle_3b_finalise, co_angle_3b_finalise, & co_distance_2b_finalise, cosnx_finalise, trihis_finalise, water_monomer_finalise, water_dimer_finalise, rdf_finalise, as_distance_2b_finalise, alex_finalise, & A2_dimer_finalise, AB_dimer_finalise, atom_real_space_finalise, power_so3_finalise, power_SO4_finalise, soap_finalise, distance_Nb_finalise, soap_turbo_finalise, & - AN_monomer_finalise, general_monomer_finalise, general_dimer_finalise, general_trimer_finalise, molecule_lo_d_finalise, com_dimer_finalise, & + AN_monomer_finalise, general_monomer_finalise, general_dimer_finalise, general_trimer_finalise, water_trimer_finalise, molecule_lo_d_finalise, com_dimer_finalise, & bond_real_space_finalise, soap_express_finalise #else module procedure descriptor_finalise, descriptor_data_finalise, RadialFunction_finalise, fourier_so4_finalise, cplx_2d_array1_finalise, cplx_3d_array2_finalise, & @@ -612,7 +614,7 @@ module descriptors_module co_distance_2b_calc, cosnx_calc, trihis_calc, water_monomer_calc, water_dimer_calc, A2_dimer_calc, AB_dimer_calc, atom_real_space_calc, & power_so3_calc, power_SO4_calc, soap_calc, rdf_calc, as_distance_2b_calc, & distance_Nb_calc, alex_calc, soap_turbo_calc, & - AN_monomer_calc, soap_express_calc, general_monomer_calc, general_dimer_calc, general_trimer_calc, molecule_lo_d_calc, com_dimer_calc, bond_real_space_calc + AN_monomer_calc, soap_express_calc, general_monomer_calc, general_dimer_calc, general_trimer_calc, water_trimer_calc, molecule_lo_d_calc, com_dimer_calc, bond_real_space_calc #else module procedure descriptor_calc, descriptor_calc_array, bispectrum_SO4_calc, bispectrum_SO3_calc, behler_calc, distance_2b_calc, coordination_calc, angle_3b_calc, co_angle_3b_calc, & co_distance_2b_calc, cosnx_calc, trihis_calc, water_monomer_calc, water_dimer_calc, A2_dimer_calc, AB_dimer_calc, atom_real_space_calc, & @@ -628,7 +630,7 @@ module descriptors_module module procedure descriptor_cutoff, bispectrum_SO4_cutoff, bispectrum_SO3_cutoff, behler_cutoff, distance_2b_cutoff, coordination_cutoff, angle_3b_cutoff, co_angle_3b_cutoff, & co_distance_2b_cutoff, cosnx_cutoff, trihis_cutoff, water_monomer_cutoff, water_dimer_cutoff, A2_dimer_cutoff, AB_dimer_cutoff, atom_real_space_cutoff, & power_so3_cutoff, power_SO4_cutoff, soap_cutoff, alex_cutoff, distance_Nb_cutoff, rdf_cutoff, as_distance_2b_cutoff, soap_turbo_cutoff, & - molecule_lo_d_cutoff, com_dimer_cutoff, soap_express_cutoff, AN_monomer_cutoff, general_monomer_cutoff, general_dimer_cutoff, general_trimer_cutoff, bond_real_space_cutoff + molecule_lo_d_cutoff, com_dimer_cutoff, soap_express_cutoff, AN_monomer_cutoff, general_monomer_cutoff, general_dimer_cutoff, general_trimer_cutoff, water_trimer_cutoff, bond_real_space_cutoff #else module procedure descriptor_cutoff, bispectrum_SO4_cutoff, bispectrum_SO3_cutoff, behler_cutoff, distance_2b_cutoff, coordination_cutoff, angle_3b_cutoff, co_angle_3b_cutoff, & co_distance_2b_cutoff, cosnx_cutoff, trihis_cutoff, water_monomer_cutoff, water_dimer_cutoff, A2_dimer_cutoff, AB_dimer_cutoff, atom_real_space_cutoff, & @@ -643,7 +645,7 @@ module descriptors_module co_distance_2b_sizes, cosnx_sizes, trihis_sizes, water_monomer_sizes, water_dimer_sizes, A2_dimer_sizes, AB_dimer_sizes, atom_real_space_sizes, & power_so3_sizes, power_SO4_sizes, soap_sizes, rdf_sizes, as_distance_2b_sizes, & alex_sizes, distance_Nb_sizes, soap_turbo_sizes, & - molecule_lo_d_sizes, com_dimer_sizes, soap_express_sizes, AN_monomer_sizes, general_monomer_sizes, general_dimer_sizes, general_trimer_sizes, bond_real_space_sizes + molecule_lo_d_sizes, com_dimer_sizes, soap_express_sizes, AN_monomer_sizes, general_monomer_sizes, general_dimer_sizes, general_trimer_sizes, water_trimer_sizes, bond_real_space_sizes #else module procedure descriptor_sizes, bispectrum_SO4_sizes, bispectrum_SO3_sizes, behler_sizes, distance_2b_sizes, coordination_sizes, angle_3b_sizes, co_angle_3b_sizes, & co_distance_2b_sizes, cosnx_sizes, trihis_sizes, water_monomer_sizes, water_dimer_sizes, A2_dimer_sizes, AB_dimer_sizes, atom_real_space_sizes, & @@ -676,9 +678,9 @@ function get_descriptor_type(args_str,error) logical :: is_bispectrum_so4, is_bispectrum_so3, is_behler, is_distance_2b, is_coordination, is_angle_3b, & is_co_angle_3b, is_co_distance_2b, is_cosnx, is_trihis, is_water_monomer, is_water_dimer, is_A2_dimer, & is_AB_dimer, is_bond_real_space, is_atom_real_space, is_power_so3, is_power_so4, is_soap, & - is_AN_monomer, is_general_monomer, is_general_dimer, is_general_trimer, is_rdf, is_as_distance_2b, & + is_AN_monomer, is_general_monomer, is_general_dimer, is_general_trimer, is_water_trimer, is_rdf, is_as_distance_2b, & is_molecule_lo_d, is_alex, is_com_dimer, is_distance_Nb, is_soap_express, is_soap_turbo - + integer n_true INIT_ERROR(error) call initialise(params) @@ -705,6 +707,7 @@ function get_descriptor_type(args_str,error) call param_register(params, 'general_monomer', 'false', is_general_monomer, help_string="Type of descriptor is general_monomer.") call param_register(params, 'general_dimer', 'false', is_general_dimer, help_string="Type of descriptor is general_dimer.") call param_register(params, 'general_trimer', 'false', is_general_trimer, help_string="Type of descriptor is general_trimer.") + call param_register(params, 'water_trimer', 'false', is_water_trimer, help_string="Type of descriptor is water_trimer.") call param_register(params, 'rdf', 'false', is_rdf, help_string="Type of descriptor is rdf.") call param_register(params, 'as_distance_2b', 'false', is_as_distance_2b, help_string="Type of descriptor is as_distance_2b.") call param_register(params, 'molecule_lo_d', 'false', is_molecule_lo_d, help_string="Type of descriptor is molecule_lo_d.") @@ -718,12 +721,13 @@ function get_descriptor_type(args_str,error) RAISE_ERROR("descriptor_initialise failed to parse args_str='"//trim(args_str)//"'", error) endif call finalise(params) - - if (count( (/is_bispectrum_so4, is_bispectrum_so3, is_behler, is_distance_2b, is_coordination, is_angle_3b, is_co_angle_3b, is_co_distance_2b, & + + n_true = count( (/is_bispectrum_so4, is_bispectrum_so3, is_behler, is_distance_2b, is_coordination, is_angle_3b, is_co_angle_3b, is_co_distance_2b, & is_cosnx, is_trihis, is_water_monomer, is_water_dimer, is_A2_dimer, is_AB_dimer, is_bond_real_space, is_atom_real_space, is_power_so3, is_power_so4, & - is_soap, is_AN_monomer, is_general_monomer, is_general_dimer, is_general_trimer, is_rdf, is_as_distance_2b, is_molecule_lo_d, is_alex, is_com_dimer, & - is_distance_Nb, is_soap_express, is_soap_turbo /) ) /= 1) then - RAISE_ERROR("descriptor_initialise found too few or too many IP Model types args_str='"//trim(args_str)//"'", error) + is_soap, is_AN_monomer, is_general_monomer, is_general_dimer, is_general_trimer, is_water_trimer, is_rdf, is_as_distance_2b, is_molecule_lo_d, is_alex, is_com_dimer, & + is_distance_Nb, is_soap_express, is_soap_turbo /) ) + if (n_true/= 1) then + RAISE_ERROR("descriptor_initialise found "//n_true//" IP Model types args_str='"//trim(args_str)//"'", error) endif get_descriptor_type = DT_NONE @@ -774,6 +778,8 @@ function get_descriptor_type(args_str,error) get_descriptor_type = DT_GENERAL_DIMER elseif( is_general_trimer ) then get_descriptor_type = DT_GENERAL_TRIMER + elseif( is_water_trimer ) then + get_descriptor_type = DT_WATER_TRIMER elseif( is_rdf ) then get_descriptor_type = DT_RDF elseif( is_as_distance_2b ) then @@ -867,6 +873,8 @@ subroutine descriptor_initialise(this,args_str,error) call initialise(this%descriptor_general_dimer,args_str,error) case(DT_GENERAL_TRIMER) call initialise(this%descriptor_general_trimer,args_str,error) + case(DT_WATER_TRIMER) + call initialise(this%descriptor_water_trimer,args_str,error) case(DT_SOAP_EXPRESS) call initialise(this%descriptor_soap_express,args_str,error) #endif @@ -938,6 +946,8 @@ subroutine descriptor_finalise(this,error) call finalise(this%descriptor_general_dimer,error) case(DT_GENERAL_TRIMER) call finalise(this%descriptor_general_trimer,error) + case(DT_WATER_TRIMER) + call finalise(this%descriptor_water_trimer,error) case(DT_SOAP_EXPRESS) call finalise(this%descriptor_soap_express,error) case(DT_SOAP_TURBO) @@ -3690,7 +3700,7 @@ subroutine descriptor_str_add_species(this,species,descriptor_str,error) enddo enddo enddo - case(DT_GENERAL_MONOMER,DT_GENERAL_DIMER,DT_GENERAL_TRIMER,DT_WATER_MONOMER,DT_WATER_DIMER,DT_A2_DIMER,DT_AB_DIMER,DT_TRIHIS,DT_BOND_REAL_SPACE,DT_ATOM_REAL_SPACE,DT_AN_MONOMER) + case(DT_GENERAL_MONOMER,DT_GENERAL_DIMER,DT_GENERAL_TRIMER,DT_WATER_TRIMER,DT_WATER_MONOMER,DT_WATER_DIMER,DT_A2_DIMER,DT_AB_DIMER,DT_TRIHIS,DT_BOND_REAL_SPACE,DT_ATOM_REAL_SPACE,DT_AN_MONOMER) allocate(descriptor_str(1)) descriptor_str(1) = trim(this) case(DT_DISTANCE_NB) @@ -3842,6 +3852,8 @@ subroutine descriptor_calc(this,at,descriptor_out,do_descriptor,do_grad_descript call calc(this%descriptor_general_dimer,at,descriptor_out,do_descriptor,do_grad_descriptor,args_str,error) case(DT_GENERAL_TRIMER) call calc(this%descriptor_general_trimer,at,descriptor_out,do_descriptor,do_grad_descriptor,args_str,error) + case(DT_WATER_TRIMER) + call calc(this%descriptor_water_trimer,at,descriptor_out,do_descriptor,do_grad_descriptor,args_str,error) case(DT_MOLECULE_LO_D) call calc(this%descriptor_molecule_lo_d,at,descriptor_out,do_descriptor,do_grad_descriptor,args_str,error) case(DT_SOAP_EXPRESS) @@ -10830,6 +10842,8 @@ function descriptor_dimensions(this,error) descriptor_dimensions = general_dimer_dimensions(this%descriptor_general_dimer,error) case(DT_GENERAL_TRIMER) descriptor_dimensions = general_trimer_dimensions(this%descriptor_general_trimer,error) + case(DT_WATER_TRIMER) + descriptor_dimensions = water_trimer_dimensions(this%descriptor_water_trimer,error) case(DT_MOLECULE_LO_D) descriptor_dimensions = molecule_lo_d_dimensions(this%descriptor_molecule_lo_d,error) case(DT_SOAP_EXPRESS) @@ -11292,6 +11306,8 @@ function descriptor_cutoff(this,error) descriptor_cutoff = cutoff(this%descriptor_general_dimer,error) case(DT_GENERAL_TRIMER) descriptor_cutoff = cutoff(this%descriptor_general_trimer,error) + case(DT_WATER_TRIMER) + descriptor_cutoff = cutoff(this%descriptor_water_trimer,error) case(DT_COM_DIMER) descriptor_cutoff = cutoff(this%descriptor_com_dimer,error) case(DT_SOAP_EXPRESS) @@ -11751,6 +11767,9 @@ subroutine descriptor_sizes(this,at,n_descriptors,n_cross,mask,n_index,error) case(DT_GENERAL_TRIMER) call general_trimer_sizes(this%descriptor_general_trimer,at, & n_descriptors,n_cross,mask=mask,n_index=n_index,error=error) + case(DT_WATER_TRIMER) + call water_trimer_sizes(this%descriptor_water_trimer,at, & + n_descriptors,n_cross,mask=mask,n_index=n_index,error=error) case(DT_COM_DIMER) call com_dimer_sizes(this%descriptor_com_dimer,at, & n_descriptors,n_cross,mask=mask,n_index=n_index,error=error) @@ -12224,7 +12243,7 @@ subroutine water_dimer_sizes(this,at,n_descriptors,n_cross,mask,n_index,error) n_descriptors = 0 n_cross = 0 -call print("mask present ? "//present(mask)) +call print("mask present ? "//present(mask),PRINT_NERD) do i = 1, at%N if(at%Z(i) == 8) then if(present(mask)) then @@ -12686,6 +12705,11 @@ function descriptor_n_permutations(this,error) RAISE_ERROR("descriptor_n_permutations: permutation_data not initialised "//this%descriptor_type,error) end if descriptor_n_permutations = this%descriptor_general_trimer%permutation_data%n_perms + case(DT_WATER_TRIMER) + if (.not. this%descriptor_water_trimer%permutation_data%initialised)then + RAISE_ERROR("descriptor_n_permutations: permutation_data not initialised "//this%descriptor_type,error) + end if + descriptor_n_permutations = this%descriptor_water_trimer%permutation_data%n_perms case(DT_MOLECULE_LO_D) if (.not. this%descriptor_molecule_lo_d%permutation_data%initialised)then RAISE_ERROR("descriptor_n_permutations: permutation_data not initialised "//this%descriptor_type,error) diff --git a/descriptors_noncommercial.inc b/descriptors_noncommercial.inc index d7cb9bbe..5d82c92a 100644 --- a/descriptors_noncommercial.inc +++ b/descriptors_noncommercial.inc @@ -1097,6 +1097,9 @@ deallocate(intermolecular) endsubroutine general_trimer_initialise + + + subroutine general_trimer_finalise(this,error) type(general_trimer), intent(inout) :: this integer, optional, intent(out) :: error @@ -1121,7 +1124,6 @@ endsubroutine general_trimer_finalise - subroutine molecule_lo_d_initialise(this,args_str,error) type(molecule_lo_d), intent(inout) :: this character(len=*), intent(in) :: args_str @@ -2571,15 +2573,15 @@ integer :: d, n_descriptors, n_cross, dimer_size, trimer_size, & i, j, k, n, m, i_atomic, j_atomic, start, finish, i_desc,this_diff,triplet_pos, & n_index - integer :: monomer_one_size, monomer_two_size, monomer_three_size, n_monomer_one, n_monomer_two, n_monomer_three, n_products, contributor_loc + integer :: monomer_one_size, monomer_two_size, monomer_three_size, n_monomer_one, n_monomer_two, n_monomer_three, n_products!, contributor_loc integer, dimension(1) :: unit_array - real(dp) :: temp_dist + real(dp) :: temp_dist, desc_der real(dp), dimension(3) :: diff_one_two, diff_one_three,diff_two_three,mean_pos_one,mean_pos_two,mean_pos_three - real(dp), dimension(3) :: grad_cut_one, grad_cut_two, grad_cut_three + real(dp), dimension(3) :: grad_cut_one, grad_cut_two, grad_cut_three, unit_vec real(dp) :: dist_one_two, dist_one_three, dist_two_three, cut12, cut13, cut23, dcut12, dcut13, dcut23 real(dp), dimension(:), allocatable :: dist_vec, primitive_cutoffs, primitive_cutoff_grads real(dp), dimension(:,:), allocatable :: interatomic_distances,diffs_one,diffs_two,diffs_three,triplets_diffs - real(dp), dimension(:,:,:), allocatable :: interatomic_vectors + real(dp), dimension(:,:,:), allocatable :: interatomic_vectors, unit_vectors integer, dimension(3) :: temp_shift, shift_one_two,shift_one_three integer, dimension(:), allocatable :: atomic_index_dimer, atomic_index_trimer, atomic_index_one, atomic_index_two, atomic_index_three,triplets_diffs_map integer, dimension(:,:), allocatable :: monomer_one_index, monomer_two_index, monomer_three_index, monomer_triplets @@ -2644,6 +2646,7 @@ allocate(atomic_index_dimer(dimer_size)) allocate(atomic_index_trimer(trimer_size)) + allocate(unit_vectors(d,3,trimer_size)) allocate(interatomic_vectors(trimer_size,trimer_size,3)) allocate(interatomic_distances(trimer_size,trimer_size)) allocate(associated_to_monomer(at%N)) @@ -2956,16 +2959,25 @@ call print("TRIM 9",PRINT_NERD) call print('TRIM Lattice="10.0000000 0.00000000 0.00000000 0.00000000 10.0000000 0.00000000 0.00000000 0.00000000 10.00000" Properties=Z:I:1:pos:R:3',PRINT_NERD) do i_atomic=1,trimer_size - call print("TRIM "//at%Z(atomic_index_trimer(i_atomic))//" "//descriptor_out%x(i_desc)%pos(:,i_atomic),PRINT_NERD) + call print("TRIM "//at%Z(atomic_index_trimer(i_atomic))//" "//descriptor_out%x(i_desc)%pos(:,i_atomic),PRINT_NERD) end do - - !build the grad_data matrix + + ! Build the grad_data matrix + unit_vectors = 0 descriptor_out%x(i_desc)%has_grad_data(:) = .true. do k=1,d i_atomic = this%component_atoms(k,1) j_atomic = this%component_atoms(k,2) - descriptor_out%x(i_desc)%grad_data(k,:,i_atomic) = - this%power * (dist_vec(k)+this%dist_shift)**(this%power-1.0_dp) * interatomic_vectors(i_atomic,j_atomic,:) / interatomic_distances(i_atomic,j_atomic) ! descriptor wrt atom i_atomic - descriptor_out%x(i_desc)%grad_data(k,:,j_atomic) = -descriptor_out%x(i_desc)%grad_data(k,:,i_atomic) ! descriptor wrt j_atomic + + unit_vec = interatomic_vectors(i_atomic,j_atomic,:) / interatomic_distances(i_atomic,j_atomic) + desc_der = this%power * (dist_vec(k)+this%dist_shift)**(this%power-1.0_dp) + + descriptor_out%x(i_desc)%grad_data(k,:,i_atomic) = - desc_der * unit_vec ! descriptor wrt atom i_atomic + descriptor_out%x(i_desc)%grad_data(k,:,j_atomic) = desc_der * unit_vec ! descriptor wrt j_atomic + + ! Unit vectors for cutoff derivatives + unit_vectors(k,:,i_atomic) = -unit_vec + unit_vectors(k,:,j_atomic) = unit_vec end do cutoff_type_grad: if (this%mpifind) then @@ -3006,8 +3018,8 @@ if (this%cutoff_contributor(m)) then if (any(this%component_atoms(k,:) == i_atomic) .or. any(this%component_atoms(m,:) == i_atomic)) then descriptor_out%x(i_desc)%grad_covariance_cutoff(:,i_atomic) = descriptor_out%x(i_desc)%grad_covariance_cutoff(:,i_atomic) & - + primitive_cutoffs(m)*primitive_cutoff_grads(k)*descriptor_out%x(i_desc)%grad_data(k,:,i_atomic) & - + primitive_cutoffs(k)*primitive_cutoff_grads(m)*descriptor_out%x(i_desc)%grad_data(m,:,i_atomic) + + primitive_cutoffs(m) * primitive_cutoff_grads(k) * unit_vectors(k,:,i_atomic) & + + primitive_cutoffs(k) * primitive_cutoff_grads(m) * unit_vectors(m,:,i_atomic) end if end if end do @@ -3018,9 +3030,9 @@ !normalisation factor descriptor_out%x(i_desc)%grad_covariance_cutoff(:,:) = descriptor_out%x(i_desc)%grad_covariance_cutoff(:,:) / n_products end if cutoff_type_grad - + end if calc_grad_descriptor - + end do loop_trimers end do loop_triplets end do loop_monomer_pairs @@ -3057,7 +3069,6 @@ endsubroutine general_trimer_calc - subroutine molecule_lo_d_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,args_str,error) type(molecule_lo_d), intent(in) :: this type(atoms), intent(in) :: at @@ -3397,6 +3408,8 @@ endfunction general_trimer_dimensions + + function molecule_lo_d_dimensions(this,error) result(i) type(molecule_lo_d), intent(in) :: this integer, optional, intent(out) :: error @@ -3493,6 +3506,7 @@ endfunction general_trimer_cutoff + function molecule_lo_d_cutoff(this,error) type(molecule_lo_d), intent(in) :: this integer, optional, intent(out) :: error @@ -3803,6 +3817,7 @@ endsubroutine general_trimer_sizes + subroutine molecule_lo_d_sizes(this,at,n_descriptors,n_cross,mask,n_index,error) type(molecule_lo_d), intent(in) :: this type(atoms), intent(in) :: at @@ -4826,4 +4841,452 @@ !!!!!!!! - +!!!!!!!! +!!! 3-body water descriptor starts here +!!! -- Jonatan Öström, @sujona, jonatan.ostrom@gmail.com +!!! +!!! This is an adaptation of the general_trimer routines to calculate water trimers faster. +!!! These routines (should) give equivalent results as long as +!!! - the lattice is orthogonal +!!! - the shortest lattice vector is > 2*trimer_cutoff +!!! - there is only water in the simulation +!!! The main issues of the general_trimer (and dimer) routines are that large arrays (monomer_triplets, triplets_diffs, ...) are +!!! - reallocated at each found trimer (or dimer) +!!! - double-checked to avoid double inclusions +!!!!!!!! + + + +subroutine water_trimer_initialise(this,args_str,error) + type(water_trimer), intent(inout) :: this + character(len=*), intent(in) :: args_str + integer, optional, intent(out) :: error + integer :: i,j, start, finish + integer, parameter :: trimer_size = 9, d = 36, at_per_mol = 3 + integer, parameter :: signature(*) = [8,1,1,8,1,1,8,1,1] + logical :: intermolecular(9,9), ordercheck + type(Dictionary) :: params + + + INIT_ERROR(error) + + call finalise(this) + + call initialise(params) + call param_register(params, 'cutoff', '4.00', this%cutoff, help_string="Cutoff(intermolecular) for water_trimer-type descriptors") + call param_register(params, 'monomer_cutoff', '1.90', this%monomer_cutoff, help_string="Monomer cutoff for water_trimer-type descriptors") + call param_register(params, 'cutoff_transition_width', '1.00', this%cutoff_transition_width, help_string="Width of smooth cutoff region for water_trimer-type descriptors") + call param_register(params, 'atom_ordercheck', 'false', ordercheck, help_string="T: find molecules. F: go by order of atoms") + call param_register(params, 'strict', 'true', this%strict, help_string="Raise error if not all atoms assigned to monomer or if no monomer pairs found") + call param_register(params, 'power', '1.0', this%power, help_string="Power of distances to be used in the kernel") + call param_register(params, 'dist_shift', '0.0', this%dist_shift, help_string="Distance shift for inverse distance descriptors.") + + + if (ordercheck)then + RAISE_ERROR("This rouine doesn't do ordercheck, OHHOHHOHH is necessary",error) + endif + if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='water_trimer_initialise args_str')) then + RAISE_ERROR("water_trimer_initialise failed to parse args_str='"//trim(args_str)//"'", error) + endif + call finalise(params) + + call permutation_data_initialise(this%permutation_data,signature_one=this%signature,signature_two=this%signature,signature_three=this%signature,internal_swaps_only=this%internal_swaps_only,error=error) + + + ! allocate(signature(trimer_size)) + ! allocate(intermolecular(trimer_size,trimer_size)) + + ! allocate(this%is_intermolecular(d)) + ! allocate(this%cutoff_contributor(d)) + ! allocate(this%component_atoms(d,2)) + + ! signature(1:n_atoms_one) = this%signature_one + ! signature(1+n_atoms_one:n_atoms_one+n_atoms_two) = this%signature_two + ! signature(1+n_atoms_one+n_atoms_two:trimer_size) = this%signature_three + intermolecular = .false. + this%cutoff_contributor=.false. + + ! Set the intermolecular interactions + do i=1,at_per_mol + do j=1+at_per_mol,trimer_size + intermolecular(i,j)=.true. + intermolecular(j,i)=.true. + end do + end do + do i=1+at_per_mol,at_per_mol+at_per_mol + do j=1+at_per_mol+at_per_mol,trimer_size + intermolecular(i,j)=.true. + intermolecular(j,i)=.true. + end do + end do + + ! Linearize intermolecular matrix, and component_atoms + start = 0 + finish=trimer_size-1 + do i=1,trimer_size + do j=1,finish-start + this%is_intermolecular(start+j) = intermolecular(i,i+j) + this%component_atoms(start+j,:) = (/ i, i+j /) + end do + start = finish + finish=finish + trimer_size-i-1 + end do + + ! Only non-hydrogens contribute to general cutoff + do i=1,d + if (this%is_intermolecular(i)) then + if (.not. signature(this%component_atoms(i,1))==1 ) then + if (.not. signature(this%component_atoms(i,2))==1 ) then + this%cutoff_contributor(i)=.true. + end if + end if + end if + end do + + this%initialised = .true. +endsubroutine water_trimer_initialise + + +subroutine water_trimer_finalise(this,error) + type(water_trimer), intent(inout) :: this + integer, optional, intent(out) :: error + INIT_ERROR(error) + if(.not. this%initialised) return + this%cutoff = 0.0_dp + this%cutoff_transition_width = 0.0_dp + this%atom_ordercheck = .false. + this%internal_swaps_only = .true. + this%use_smooth_cutoff = .false. + this%power = 1.0_dp + this%dist_shift = 0.0_dp + this%initialised = .false. +endsubroutine water_trimer_finalise + + + +function water_trimer_dimensions(this,error) result(i) + type(water_trimer), intent(in) :: this + integer, optional, intent(out) :: error + integer :: i + + INIT_ERROR(error) + + if(.not. this%initialised) then + RAISE_ERROR("water_trimer_dimensions: descriptor object not initialised", error) + endif + if(.not. this%permutation_data%initialised) then + RAISE_ERROR("water_trimer_dimensions: descriptor object's permutation data not initialised", error) + endif + + i = size(this%permutation_data%dist_vec) + +endfunction water_trimer_dimensions + +function water_trimer_cutoff(this,error) + type(water_trimer), intent(in) :: this + integer, optional, intent(out) :: error + real(dp) :: water_trimer_cutoff + + INIT_ERROR(error) + + if(.not. this%initialised) then + RAISE_ERROR("water_trimer_cutoff: descriptor object not initialised", error) + endif + + water_trimer_cutoff = this%cutoff + +endfunction water_trimer_cutoff + + + +subroutine water_trimer_sizes(this,at,n_descriptors,n_cross,mask,n_index,error) + use find_water_triplets, only: find_triplets_lars + type(water_trimer), intent(in) :: this + type(atoms), intent(in) :: at + integer, intent(out) :: n_descriptors, n_cross + logical, dimension(:), intent(in), optional :: mask + integer, intent(out), optional :: n_index + integer, optional, intent(out) :: error + + integer, dimension(:,:), allocatable :: monomer_one_index, monomer_triplets + real(dp), dimension(:,:), allocatable :: triplets_diffs + real(dp) box(3) + + INIT_ERROR(error) + + if(.not. this%initialised) then + RAISE_ERROR("water_trimer_sizes: descriptor object not initialised", error) + endif + + call water_trimer_system_check(this,at,box,error) + call find_triplets_lars(at%pos,at%N/3,box,this%cutoff,monomer_triplets,triplets_diffs,n_descriptors) + + n_cross = n_descriptors * 9 + + if(allocated(monomer_triplets)) deallocate(monomer_triplets) + if(allocated(triplets_diffs)) deallocate(triplets_diffs) + + if( present(n_index) ) n_index = 9 + +endsubroutine water_trimer_sizes + +subroutine water_trimer_system_check(this,at,box,error) + type(water_trimer), intent(in) :: this + type(atoms), intent(in) :: at + real(dp), intent(out) :: box(3) + integer, optional, intent(inout) :: error + real(dp) :: lattice(3,3) + logical :: ok_ohh, ok_box, ok_ort + integer i + + lattice = at%lattice + + ! Check OHH + ok_ohh = all(at%Z(1::3) == 8) .and. all(at%Z(2::3) == 1) .and. all(at%Z(3::3) == 1) + + do i = 1,3 + box(i) = lattice(i,i) + lattice(i,i) = 0 + enddo + + ! Check orthogonal + ok_ort = all(abs(lattice)<1d-14) .and. at%is_orthorhombic + + ! Fix box (lattice) size if it is 0 + ! If no Lattice is given: in quip (gap_fit) the box is set to something random but big enough, in quippy to 0. + if (all(box<1)) box = 1000 ! If box is 0 (<1) set it to very big + + + ! Check that box is large enough + ok_box = all(box > 2*this%cutoff) + + if (.not.all([ok_ohh,ok_ort,ok_box]))then + ! STOP __FILE__//":"//__LINE__//": water_trimer_calc: Inapplicable: ohh_order="//ok_ohh//" orthogonal="//ok_ort//" large_enough="//ok_box//" this box="//box + RAISE_ERROR("water_trimer: Inapplicable: ohh_order="//ok_ohh//" orthogonal="//ok_ort//" large_enough="//ok_box, error) + endif +end subroutine + + +subroutine water_trimer_calc(this,at,descriptor_out,do_descriptor,do_grad_descriptor,args_str,error) !////////////////////////////////////////////////////////////////////////////////// + use find_water_triplets, only: find_triplets_lars + type(water_trimer), intent(in) :: this + type(atoms), intent(in) :: at + type(descriptor_data), intent(out) :: descriptor_out + logical, intent(in), optional :: do_descriptor, do_grad_descriptor + character(len=*), intent(in), optional :: args_str + integer, optional, intent(out) :: error + integer :: n_descriptors, i, i_desc + real(dp), dimension(:,:), allocatable :: triplets_diffs + integer, dimension(:,:), allocatable :: itriplets + logical :: my_do_descriptor, my_do_grad_descriptor + integer, parameter :: trimer_size = 9, nd = 36 + integer, dimension(3,at%N/3) :: imonomers + real(dp) :: box(3) + + INIT_ERROR(error) + + call system_timer('water_trimer_calc') + + if(.not. this%initialised) then + RAISE_ERROR("water_trimer_calc: descriptor object not initialised", error) + endif + + my_do_descriptor = optional_default(.false., do_descriptor) + my_do_grad_descriptor = .true. !optional_default(.false., do_grad_descriptor) ! fix!!! + + if( .not. my_do_descriptor .and. .not. my_do_grad_descriptor ) return + + call finalise(descriptor_out) + + + call water_trimer_system_check(this,at,box,error) + call find_triplets_lars(at%pos,at%N/3,box,this%cutoff,itriplets,triplets_diffs,n_descriptors) + + imonomers = reshape([(i,i=1,at%N)],[3,at%N/3]) + + ! Print stuff? + call system_timer('water_trimer_calc: find_monomer_triplets') + call print("itriplets("//size(itriplets,1)//", "//size(itriplets,2)//")", PRINT_NERD) + call print("triplets_diffs("//size(triplets_diffs,1)//", "//size(triplets_diffs,2)//")", PRINT_NERD) + call print("ready to construct "//n_descriptors //" descriptors",PRINT_NERD) + + ! Allocate + allocate(descriptor_out%x(n_descriptors)) + do i = 1, n_descriptors + if(my_do_descriptor) then + allocate(descriptor_out%x(i)%data(nd)) + allocate(descriptor_out%x(i)%ci(trimer_size)) + descriptor_out%x(i)%has_data = .false. + descriptor_out%x(i)%data = 0.0_dp + descriptor_out%x(i)%ci = 0 + descriptor_out%x(i)%covariance_cutoff = 1.0_dp + endif + if(my_do_grad_descriptor) then + allocate(descriptor_out%x(i)%grad_data(nd,3,trimer_size)) + allocate(descriptor_out%x(i)%ii(trimer_size)) + allocate(descriptor_out%x(i)%pos(3,trimer_size)) + allocate(descriptor_out%x(i)%has_grad_data(trimer_size)) + allocate(descriptor_out%x(i)%grad_covariance_cutoff(3,trimer_size)) + descriptor_out%x(i)%has_grad_data = .false. + descriptor_out%x(i)%grad_data = 0.0_dp + descriptor_out%x(i)%ii = 0 + descriptor_out%x(i)%pos = 0.0_dp + descriptor_out%x(i)%grad_covariance_cutoff = 0.0_dp + endif + enddo + + ! Loop through trimers + ! Could do OpenMP here but that's done externally for the training... Hmm but SOAP has it right? + ! Blocks are in place, so making functions or openMP is easy. + the_loop : do i_desc = 1, n_descriptors + + block + integer, parameter :: monomer_size = 3, trimer_size = 9, nd = 36 + real(dp) cut + real(dp), dimension(3,9) :: pos, gcuts + real(dp), dimension(nd) :: desc0, rlen + ! real(dp), dimension(3,36) :: rvec + ! real(dp), dimension(3,9,36) :: gdesc0 + real(dp), dimension(nd,3,9) :: gdesc0 + integer ii, jj + integer i1, i2,i3 + real(dp) :: diffs1(3,monomer_size) + real(dp) :: diffs2(3,monomer_size) + real(dp) :: diffs3(3,monomer_size) + integer :: iatoms1(monomer_size) + integer :: iatmos2(monomer_size) + integer :: iatoms3(monomer_size) + integer :: iatoms123(trimer_size) + real(dp), dimension(3) :: diff12, diff13 + real(dp), dimension(3) :: mean_pos1,mean_pos2,mean_pos3 + + + i1 = itriplets(1,i_desc) + i2 = itriplets(2,i_desc) + i3 = itriplets(3,i_desc) + + iatoms1 = imonomers(:,i1) !store the indices of atoms in this monomer + iatmos2 = imonomers(:,i2) + iatoms3 = imonomers(:,i3) + + iatoms123=(/iatoms1,iatmos2,iatoms3/) + + mean_pos1 = calc_mean_pos(at,iatoms1) + mean_pos2 = calc_mean_pos(at,iatmos2) + mean_pos3 = calc_mean_pos(at,iatoms3) + + + diff12 = triplets_diffs(1:3,i_desc) + diff13 = triplets_diffs(4:6,i_desc) + + + do ii=1,3 + diffs1(:,ii) = diff_min_image(at,at%pos(:,iatoms1(ii)),mean_pos1) + diffs2(:,ii) = diff_min_image(at,at%pos(:,iatmos2(ii)),mean_pos2) + diffs3(:,ii) = diff_min_image(at,at%pos(:,iatoms3(ii)),mean_pos3) + end do + + + do ii=1,3 + pos(:, ii) = mean_pos1 - diffs1(:,ii) + pos(:,3+ii) = mean_pos1 + diff12 - diffs2(:,ii) + pos(:,6+ii) = mean_pos1 + diff13 - diffs3(:,ii) + end do + + ! Calcualate descriptor + block + real(dp) a0, shft, pow + real(dp), dimension(3) :: a3, b3 + integer ii, jj, ij + + shft = this%dist_shift + pow = this%power + gdesc0 = 0 + do ij = 1, nd + jj = this%component_atoms(ij,1) + ii = this%component_atoms(ij,2) + a3 = pos(:,ii)-pos(:,jj) + a0 = norm(a3) + ! rvec(:,ij) = a3 + rlen(ij) = a0 + desc0(ij) = (a0 + shft) ** pow + + ! Gradient + b3 = - pow * (a0+shft)**(pow-1.0_dp) * a3 / a0 + gdesc0(ij,:,jj) = b3 + gdesc0(ij,:,ii) = -b3 + + enddo + endblock + + ! Calculate cutoff + block + real(dp), dimension(3) :: oo12, oo13, oo23 + real(dp), dimension(3) :: gcut1, gcut2, gcut3 + real(dp) :: o12, o13, o23 + real(dp) :: cut12, cut13, cut23 + real(dp) :: dcut12, dcut13, dcut23 + real(dp) :: cutoff, width + cutoff= this%cutoff + width = this%cutoff_transition_width + + oo12 = pos(:,4)-pos(:,1) + oo13 = pos(:,7)-pos(:,1) + oo23 = pos(:,7)-pos(:,4) + o12 = norm(oo12) + o13 = norm(oo13) + o23 = norm(oo23) + cut12 = coordination_function(o12,cutoff,width) + cut13 = coordination_function(o13,cutoff,width) + cut23 = coordination_function(o23,cutoff,width) + cut = (cut12*cut13 + cut12*cut23 + cut13*cut23) / 3 + + ! Gradient + dcut12 = dcoordination_function(o12, cutoff, width) + dcut13 = dcoordination_function(o13, cutoff, width) + dcut23 = dcoordination_function(o23, cutoff, width) + gcut1 = (- oo12 / o12 * dcut12 * (cut13+cut23) - oo13 / o13 * dcut13 * (cut12+cut23)) / 3 + gcut2 = ( oo12 / o12 * dcut12 * (cut13+cut23) - oo23 / o23 * dcut23 * (cut12+cut13)) / 3 + gcut3 = ( oo13 / o13 * dcut13 * (cut12+cut23) + oo23 / o23 * dcut23 * (cut12+cut13)) / 3 + + gcuts = 0 + gcuts(:,1) = gcut1 + gcuts(:,4) = gcut2 + gcuts(:,7) = gcut3 + endblock + + if(my_do_descriptor) then + descriptor_out%x(i_desc)%has_data = .true. + descriptor_out%x(i_desc)%data(:) = desc0 + descriptor_out%x(i_desc)%ci(:) = iatoms123 + descriptor_out%x(i_desc)%covariance_cutoff = cut + end if + + if(my_do_grad_descriptor) then + descriptor_out%x(i_desc)%has_grad_data = .true. + descriptor_out%x(i_desc)%grad_data(:,:,:) = gdesc0 + descriptor_out%x(i_desc)%ii(:) = iatoms123 + descriptor_out%x(i_desc)%pos(:,:) = pos + descriptor_out%x(i_desc)%grad_covariance_cutoff(:,:) = gcuts + end if + + endblock + if (.false.) then ! to compare with general trimer + write(12,'(a,*(f10.5))') " desc = ", descriptor_out%x(i_desc)%data + write(12,'(a,*(f10.5))') "g desc = ", descriptor_out%x(i_desc)%grad_data + write(12,'(a,*(f10.5))') " cut = ", descriptor_out%x(i_desc)%covariance_cutoff + write(12,'(a,*(f10.5))') "g cut = ", descriptor_out%x(i_desc)%grad_covariance_cutoff + stop "testing" + endif + enddo the_loop + + if (allocated(itriplets)) deallocate(itriplets) + if (allocated(triplets_diffs)) deallocate(triplets_diffs) + + call system_timer('water_trimer_calc') + +endsubroutine water_trimer_calc !//////////////////////////////////////////////////////////////////////////////////////////////////////////// + +!!!!!!!! +!!! 3-body water descriptor ends here +!!!!!!!! diff --git a/descriptors_noncommercial_permutations.inc b/descriptors_noncommercial_permutations.inc index 55c6b805..28a15e99 100644 --- a/descriptors_noncommercial_permutations.inc +++ b/descriptors_noncommercial_permutations.inc @@ -126,6 +126,21 @@ permutations=my_permutation_data%dist_vec_permutations + case(DT_WATER_TRIMER) + if (.not. this%descriptor_water_trimer%permutation_data%initialised)then + RAISE_ERROR("descriptor_permutations: permutation_data not initialised "//this%descriptor_type,error) + else if (this%descriptor_water_trimer%permutation_data%perm_number /= 1) then + RAISE_ERROR("descriptor_permutations: permutation_data%perm_number must be initialised to one"//this%descriptor_type,error) + end if + + call permutation_data_copy(my_permutation_data, this%descriptor_water_trimer%permutation_data) + + if (my_permutation_data%n_perms > 1) then + call next(my_permutation_data, 1) + end if + + permutations=my_permutation_data%dist_vec_permutations + case(DT_MOLECULE_LO_D) if (.not. this%descriptor_molecule_lo_d%permutation_data%initialised) then RAISE_ERROR("descriptor_permutations: permutation_data not initialised "//this%descriptor_type,error) diff --git a/descriptors_noncommercial_types.inc b/descriptors_noncommercial_types.inc index 2a6d34bf..03648f9e 100644 --- a/descriptors_noncommercial_types.inc +++ b/descriptors_noncommercial_types.inc @@ -100,6 +100,20 @@ endtype general_trimer + type water_trimer + !!! 3-body water descriptor + !!! -- Jonatan Öström, @sujona, jonatan.ostrom@gmail.com + type(permutation_data_type) :: permutation_data + integer, dimension(3) :: signature = [8,1,1] + real(dp) :: cutoff, cutoff_transition_width, monomer_cutoff + logical :: atom_ordercheck = .false., internal_swaps_only = .true., use_smooth_cutoff + logical :: initialised = .false., strict + real(dp) :: power,dist_shift + integer, dimension(36,2):: component_atoms + logical, dimension(36) :: is_intermolecular, cutoff_contributor + endtype water_trimer + + type molecule_lo_d type(permutation_data_type) :: permutation_data type(Atoms) :: template_atoms diff --git a/find_water_triplets_noncommercial.f95 b/find_water_triplets_noncommercial.f95 new file mode 100644 index 00000000..e7731e6a --- /dev/null +++ b/find_water_triplets_noncommercial.f95 @@ -0,0 +1,485 @@ +! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX +! HND X +! HND X GAP (Gaussian Approximation Potental) +! HND X +! HND X +! HND X Portions of GAP were written by Albert Bartok-Partay, Gabor Csanyi, +! HND X Copyright 2006-2021. +! HND X +! HND X Portions of GAP were written by Noam Bernstein as part of +! HND X his employment for the U.S. Government, and are not subject +! HND X to copyright in the USA. +! HND X +! HND X GAP is published and distributed under the +! HND X Academic Software License v1.0 (ASL) +! HND X +! HND X GAP is distributed in the hope that it will be useful for non-commercial +! HND X academic research, but WITHOUT ANY WARRANTY; without even the implied +! HND X warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! HND X ASL for more details. +! HND X +! HND X You should have received a copy of the ASL along with this program +! HND X (e.g. in a LICENSE.md file); if not, you can write to the original licensors, +! HND X Gabor Csanyi or Albert Bartok-Partay. The ASL is also published at +! HND X http://github.com/gabor1/ASL +! HND X +! HND X When using this software, please cite the following reference: +! HND X +! HND X A. P. Bartok et al Physical Review Letters vol 104 p136403 (2010) +! HND X +! HND X When using the SOAP kernel or its variants, please additionally cite: +! HND X +! HND X A. P. Bartok et al Physical Review B vol 87 p184115 (2013) +! HND X +! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX + +!!!!!!!! +!!! This file was written by Jonatan Öström (@sujona, jonatan.ostrom@gmail.com) and Lars G.M. Pettersson, Stockholm University +!!! Here are implementations of water-trimer search routines, that work only for orthogonal unit cells with the shortest dimension longer than 2 x (3-body cutoff) +!!!!!!!! + +module find_water_triplets + ! use backend + implicit none + integer, parameter :: dp = kind(0d0) + + interface insert + module procedure insert_i2, insert_r2, insert_i1 + end interface + + +contains + +! //////////////////////////////////////////////// +! Dynamic insert in allocatable array +! insert(a,i,x) does a(i) = x but reallocates a to length 2*i if len(a) < i + +subroutine insert_i1(array,ii,val) + integer, intent(inout), allocatable :: array(:) + integer, intent(in) :: ii, val + integer, allocatable :: tmp(:) + if (ii>size(array))then + tmp = array + deallocate(array) + allocate(array(2*ii)) + array(:size(tmp)) = tmp + endif + array(ii) = val +end +subroutine insert_i2(array,ii,val) + integer, intent(inout), allocatable :: array(:,:) + integer, intent(in) :: ii, val(:) + integer, allocatable :: tmp(:,:) + if (ii>size(array,2))then + tmp = array + deallocate(array) + allocate(array(size(tmp,1),2*ii)) + array(:,:size(tmp,2)) = tmp + endif + array(:,ii) = val +end +subroutine insert_r2(array,ii,val) + real(dp), intent(inout), allocatable :: array(:,:) + real(dp), intent(in) :: val(:) + integer, intent(in) :: ii + real(dp), allocatable :: tmp(:,:) + if (ii>size(array,2))then + tmp = array + deallocate(array) + allocate(array(size(tmp,1),2*ii)) + array(:,:size(tmp,2)) = tmp + endif + array(:,ii) = val +end + +! //////////////////////////////////////////////// +! Utilitity routines + +subroutine min_img(XO,NO,i1,i2,box,lsq,x12,s12) + real(dp), intent(in) :: XO(3,NO), box(3) + integer, intent(in) :: NO, i1, i2 + real(dp), intent(out) :: lsq + real(dp), intent(out) :: x12(3) + integer, intent(out) :: s12(3) + x12 = XO(:,i1)-XO(:,i2) + s12 = nint(x12/box) + x12 = x12 - s12*box + lsq = sum(x12**2) +end + +subroutine min_img_c(XO,XC,NO,i1,i2,box,lsq,x12,c12,s12) + real(dp), intent(in) :: XO(3,NO), XC(3,NO), box(3) + integer, intent(in) :: NO, i1, i2 + real(dp), intent(out) :: lsq + real(dp), intent(out) :: x12(3), c12(3) + integer, intent(out) :: s12(3) + c12 = XC(:,i2)-XC(:,i1) + s12 = -nint(c12/box) + x12 = XO(:,i2)-XO(:,i1) + x12 = x12 + s12*box + c12 = c12 + s12*box + lsq = sum(x12**2) +end + +subroutine average_position(NW,XW,XC,XO) + integer, intent(in) :: NW + real(dp), intent(in) :: XW(3,NW*3) + real(dp), intent(out) :: XC(3,NW),XO(3,NW) + integer ii,jj,kk + do ii = 1,NW + jj = 3*(ii-1) !+1,2,3 for O,H,H + XO(:,ii) = XW(:,jj+1) + XC(:,ii) = 0 + do kk = 1,3 + XC(:,ii) = XC(:,ii) + XW(:,jj+kk)/3 + enddo + enddo +end + +function wall_time() result(tt) + integer count,count_rate + real(dp) tt + call system_clock(count,count_rate) + tt = dble(count)/count_rate +end +subroutine take(tt,text) + real(dp), intent(inout) :: tt + real(dp) :: old + character(*), intent(in) :: text + old = tt + tt = wall_time() + print'(f6.3,a)',tt-old,"s "//text +end + + +! //////////////////////////////////////////////// +! Finding triplets + +subroutine find_triplets_brute_force(XW,NW,box,rcut,n_trip) + ! find water triplets with brute force + integer, intent(in) :: NW + real(dp), intent(in) :: XW(3,NW*3), box(3), rcut + integer, intent(out) :: n_trip + real(dp) :: rcut2, d2ij, d2ik,d2jk, XC(3,NW),XO(3,NW) + integer ii,jj,kk, n_short !n_2, n_3, + real(dp),dimension(3) :: xij,xik,xjk ! Oxygen diff + real(dp),dimension(3) :: cij,cik,cjk ! Center diff + integer, dimension(3) :: sij,sik,sjk ! Shifts + rcut2 = rcut**2 + + call average_position(NW,XW,XC,XO) + + !!! $omp parallel do private(d2ij, d2ik,d2jk,n_short) reduction(+:n_2,n_3) + write(13,'(a)') ' ' + write(13,'(a)') ' TRIPLETS' + write(13,'(a)') ' ' + do ii = 1,NW + do jj = ii+1,NW + call min_img_c(XO,XC,NW,ii,jj,box,d2ij,xij,cij,sij) + do kk = jj+1,NW + call min_img_c(XO,XC,NW,ii,kk,box,d2ik,xik,cik,sik) + call min_img_c(XO,XC,NW,jj,kk,box,d2jk,xjk,cjk,sjk) + n_short = count([d2ij,d2ik,d2jk]1)then + n_trip = n_trip+1 + write(13,'(3i5,2(3i3,2x),2(3f10.5,2x))') ii,jj,kk, sij, sik, cij,cik + endif + enddo + enddo + enddo +end + +subroutine find_pairs_jona(XW,NO,box,rcut,n_pair,id2, sh2, dx2) + integer, intent(in) :: NO + real(dp), intent(in) :: XW(3,NO*3), box(3), rcut + integer, intent(out) :: n_pair + real(dp) :: rcut2, d2ij, xij(3), cij(3), XO(3,NO*3), XC(3,NO*3) + integer ii,jj, sij(3) + integer, allocatable, dimension(:,:), intent(out) :: id2, sh2 + real(dp), allocatable, dimension(:,:), intent(out) :: dx2 + allocate(id2(2,0), sh2(3,0), dx2(3,0)) + call average_position(NO,XW,XC,XO) + rcut2 = rcut**2 + n_pair = 0 + do ii = 1,NO + do jj = ii+1,NO + call min_img_c(XO,XC,NO,ii,jj,box,d2ij,xij,cij,sij) + if (d2ijii in G(ii) + do jl = 1,ncount(ii) + jj = nindex(i0+jl) + j0 = offset(jj) + sij = sh2(:,i0+jl) + cij = dx2(:,i0+jl) + if (iijj in G(ii) to get ii=a, jj=b, kk=c + do kl = jl+1,ncount(ii) + kk = nindex(i0+kl) + n_trip = n_trip + 1 + sik = sh2(:,i0+kl) + cik = dx2(:,i0+kl) + call insert(id3, n_trip, [ii,jj,kk]) + call insert(sh3, n_trip, [sij,sik]) + call insert(dx3, n_trip, [cij,cik]) + enddo + ! 2. a-b-c & 3. a-c-b: take kk>ii in G(jj)\G(ii) so that ii=a, jj=b/c, kk=c/b + do kl = 1,ncount(jj) + kk = nindex(j0+kl) + if (ii N(i) means j runs over the (N)eighbors of i + ! case 1. c-a-b-(c-) : let i N(i) and k -> N(i) + n_trip = 0 + do ii = 1, NO - 1 + i0 = ioff(ii) + i1 = ioff(ii+1) + do jl = 1, num(ii) - 1 + jj = neighbor(i0 + jl) + ! sij = sh2(:,i0 + jl) + cij = dx2(:,i0 + jl) + do kl = jl + 1, num(ii) + kk = neighbor(i0 + kl) + n_trip = n_trip + 1 + ! sik = sh2(:,i0 + kl) + cik = dx2(:,i0 + kl) + + call insert(id3, n_trip, [ii,jj,kk]) + ! call insert(sh3, n_trip, [sij,sik]) + call insert(dx3, n_trip, [cij,cik]) + + enddo + enddo + ! case 2. a-b-c : let iN(i) and k -> N(j)\N(i) + do jl = 1, num(ii) + jj = neighbor(i0 + jl) + j0 = ioff(jj) + ! sij = sh2(:,i0 + jl) + cij = dx2(:,i0 + jl) + do kl = 1, num(jj) !(1) + kk = neighbor(j0 + kl) !(2) + if (any(neighbor(i0+1:i1)==kk)) cycle + n_trip = n_trip + 1 + + ! sik = sh2(:,i0 + jl) + sh2(:,j0 + kl) + cik = dx2(:,i0 + jl) + dx2(:,j0 + kl) + + call insert(id3, n_trip, [ii,jj,kk]) + ! call insert(sh3, n_trip, [sij,sik]) + call insert(dx3, n_trip, [cij,cik]) + enddo + ! case 3. a-c-b : let i=a < k=b < j=c, so for j->N(i) find k->N(j)\N(i) + ! since k[i+1,j-1], exclude k->N(i) and include k when j->N(k) + do kk = ii + 1, jj - 1 + k0 = ioff(kk) + do jn = 1, num(kk) + if (jj /= neighbor(k0+jn)) cycle ! include j->N(k) + if (any(kk==neighbor(i0+1:i1))) cycle ! exclude k->N(i) + n_trip = n_trip + 1 + + ! sik = sh2(:,i0 + jl) - sh2(:,k0 + jn) + cik = dx2(:,i0 + jl) - dx2(:,k0 + jn) + + ! insert in i