anchor.diff:  
     integer :: nat
     integer :: i, j
  -  double precision, dimension(3) :: dipole_origin
  -  call center_of_charge(r_mol, z_mol, dipole_origin)
     nat = size(r_mol, 2)
     do i = 1, cdim
       if (constr >= 0 .and. i.eq.1) then
  
positive.code:   allocate(Aaux(ndim, nmep))
    call auxvec(r_mol, r_mep, v_mep, order, Aaux, B)
    A = matmul(Aaux, transpose(Aaux))
  
    deallocate(Aaux)
  
  end subroutine geom_mat
  
  
  
  subroutine center_of_charge(r, z, c)
    implicit none
    double precision, dimension(:,:), intent(in) :: r
    double precision, dimension(:), intent(in) :: z
    double precision, dimension(3), intent(out) :: c
  
    integer :: i, nat, zsum, zval
    c = 0.0d0
    zsum = 0
  
    nat = size(r,2)
    do i = 1, nat
        zsum = zsum + z(i)
        c(:) = c(:) + r(1:3,i) * z(i)
    enddo
    c = c / zsum
  end subroutine
  
  subroutine constr_mat(r_mol, z_mol, constr, constr_q, constr_d, ndim, cdim, A, B)
    implicit none
    double precision, intent(in), dimension(:,:) :: r_mol
    double precision, intent(in), dimension(:) :: z_mol
  
negative.code:     endif
      if (constr >= 1 .and. i.gt.1) then
        do j = 1, nat
          A(j, ndim + i) = center_of_mass(i-1) - r_mol(i-1,j)
          A(ndim + i, j) = A(j, ndim +i)
        enddo
        B(ndim + i) = constr_d(i-1)
  