Bug fixed in the determination of the spin forbidden reactions

src/Fragment.f90

@@ -61,7 +61,7 @@ module Fragment_
 	use RandomUtils_
 	use Molecule_
 	use StringList_
-	use RealList_
+	use RealVector_
 
 	use GOptionsM3C_
 
@@ -929,7 +929,7 @@ end subroutine showLnWComponents
 	!!
 	function spinAvailable( this ) result( output )
 		class(Fragment), intent(in) :: this
-		type(RealList) :: output
+		type(RealVector) :: output
 
 		real(8) :: S, Si, Sj
 		integer :: i, j

src/FragmentsListBase.f90

@@ -1,6 +1,12 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!                                                                                   !!
 !! This file is part of M3C project                                                  !!
+!!                                                                                   !!
+!! Copyright (c) 2019-2020 by authors                                                !!
+!! Authors:                                                                          !!
+!!                         * Néstor F. Aguirre (2019-2020)                           !!
+!!                           [email protected]                                    !!
+!!                                                                                   !!
 !! Copyright (c) 2013-2016 Departamento de Química                                   !!
 !!                         Universidad Autónoma de Madrid                            !!
 !!                         All rights reserved.                                      !!
@@ -66,7 +72,7 @@ module FragmentsListBase_
 	use Atom_
 	use Molecule_
 	use BlocksIFileParser_
-	use RealList_
+	use RealVector_
 
 	use GOptionsM3C_
 	use Fragment_
@@ -75,6 +81,9 @@ module FragmentsListBase_
 	implicit none
 	private
 
+	public :: &
+		spinListCoupling
+
 	type, abstract, public :: FragmentsListBase
 		type(Fragment), allocatable :: clusters(:)    !< Fragment list
 		integer, allocatable :: idSorted(:)            !< Position of the molecule sorted by mass. It is calculated in updateLabel procedure
@@ -2257,43 +2266,97 @@ function spinRange( this ) result( output )
 	end function spinRange
 
 	!>
-	!! @brief
+	!! @brief Internal. Used only by spinAvailable function
 	!!
-	function spinAvailable( this ) result( output )
-		class(FragmentsListBase), intent(in) :: this
-		type(RealList) :: output
+	function spinSpinCoupling( s1, s2 ) result( st )
+		real(8), intent(in) :: s1, s2
+		type(RealVector) :: st
 
-		real(8) :: S, Si, Sj
-		integer :: i, j
+		real(8) :: s0
+		integer :: i
+		integer :: n
 
-		call output.init()
+		n = s1+s2-abs(s1-s2)+1
+		s0 = abs(s1-s2)
 
-		if( this.nMolecules() == 1 ) then
-			S = (this.clusters(1).multiplicity-1.0_8)/2.0_8
-			call output.append( S )
-		else
-			do i=1,this.nMolecules()-1
-				Si = (this.clusters(i).multiplicity-1.0_8)/2.0_8
+		call st.init( n )
+		do i=1,n
+			call st.set( i, s0+i-1 )
+		end do
 
-				if( Si < 0.0_8 ) Si=0.0_8
+	end function spinSpinCoupling
+
+	!>
+	!! @brief Internal. Used only by spinAvailable function
+	!!
+	function spinListSpinCoupling( sl1, s2 ) result( st )
+		type(RealVector), intent(in) :: sl1
+		real(8), intent(in) :: s2
+		type(RealVector) :: st
+
+		real(8) :: s0
+		integer :: i
+		integer :: n
+
+		n = maxval(sl1.data(1:sl1.size())+s2)-minval(abs(sl1.data(1:sl1.size())-s2))+1
+		s0 = minval(abs(sl1.data(1:sl1.size())-s2))
+
+		call st.init( n, 0.0_8 )
+		do i=1,n
+			call st.set( i, s0+i-1 )
+		end do
 
-				do j=i+1,this.nMolecules()
-					Sj = (this.clusters(j).multiplicity-1.0_8)/2.0_8
-
-					if( Sj < 0.0_8 ) Sj=0.0_8
-
-					S = abs(Si-Sj)
-					do while( int(2.0*S) <= int(2.0*(Si+Sj)) )
-						call output.append( S )
-
-						S = S + 1.0_8
-					end do
-				end do
+	end function spinListSpinCoupling
+
+	!>
+	!! @brief Internal. Used only by spinAvailable function
+	!!
+	function spinListCoupling( sl ) result( st )
+		type(RealVector), intent(in) :: sl
+		type(RealVector) :: st
+
+		integer :: i
+		type(RealVector) :: sij
+
+		call st.init()
+
+		if( sl.size() == 1 ) then
+			st = sl
+		else
+			sij = spinSpinCoupling( sl.at(1), sl.at(2) )
+			st = sij
+
+			do i=3,sl.size()
+				sij = spinListSpinCoupling( sij, sl.at(i) )
+				st = sij
 			end do
 		end if
+
+	end function spinListCoupling
+
+	!>
+	!! @brief
+	!!
+	function spinAvailable( this ) result( output )
+		class(FragmentsListBase), intent(in) :: this
+		type(RealVector) :: output
+
+		integer :: i
+		real(8) :: S
+		type(RealVector) :: sl
+
+		call sl.init( this.nMolecules(), 0.0_8 )
+		do i=1,this.nMolecules()
+			S = (this.clusters(i).multiplicity-1.0_8)/2.0_8
+			if( S < 0.0_8 ) S=0.0_8
+
+			call sl.append( S )
+		end do
+
+		output = spinListCoupling( sl )
 
 		if( output.size() == 0 ) then
-			write(*,*) "### ERROR ### FragmentsListBase.spinAvailable.size() == 0", this.nMolecules(), S, Si, Sj
+			write(*,*) "### ERROR ### FragmentsListBase.spinAvailable.size() == 0", this.nMolecules(), S
 		end if
 	end function spinAvailable
 

src/Reactor.f90

@@ -63,6 +63,7 @@ module Reactor_
 	use BlocksIFileParser_
 	use RealHistogram_
 	use RealList_
+	use RealVector_
 	use IntegerVector_
 	use StringIntegerMap_
 	use StringRealMap_
@@ -71,6 +72,7 @@ module Reactor_
 
 	use GOptionsM3C_
 	use Fragment_
+	use FragmentsListBase_
 	use FragmentsList_
 	use FragmentsDB_
 
@@ -127,7 +129,7 @@ module Reactor_
 	integer, private :: internalCharge = 0
 	integer, private :: internalNTrials = 0
 ! 	real(8), private :: internalReactivesSpinRange(2) = 0.0_8
-	type(RealList), private :: internalReactivesSpinAvail
+	type(RealVector), private :: internalReactivesSpinAvail
 
 	contains
 
@@ -348,56 +350,39 @@ function isSpinForbidden( multisetPositions, current ) result( output )
 		logical :: output
 
 		real(8) :: S, Si, Sj
-		type(RealList) :: spinAvail
+		type(RealVector) :: spinAvail
 		integer :: i, j
 		class(RealListIterator), pointer :: it1, it2
 
-		call spinAvail.init()
+		call spinAvail.init( current, 0.0_8 )
 
 		if( current == 1 ) then
 			S = (FragmentsDB_instance.clusters( multisetPositions(1) ).multiplicity-1.0_8)/2.0_8
-			call spinAvail.append( S )
+			call spinAvail.set( 1, S )
 		else
-			do i=1,current-1
+			do i=1,current
 				Si = (FragmentsDB_instance.clusters( multisetPositions(i) ).multiplicity-1.0_8)/2.0_8
-
 				if( Si < 0.0_8 ) Si=0.0_8
 
-				do j=i+1,current
-					Sj = (FragmentsDB_instance.clusters( multisetPositions(j) ).multiplicity-1.0_8)/2.0_8
-
-					if( Sj < 0.0_8 ) Sj=0.0_8
-
-					S = abs(Si-Sj)
-					do while( int(2.0*S) <= int(2.0*(Si+Sj)) )
-						call spinAvail.append( S )
-						S = S + 1.0_8
-					end do
-				end do
+				call spinAvail.set( i, Si )
 			end do
+
+			spinAvail = spinListCoupling( spinAvail )
+
 		end if
 
 		if( spinAvail.size() == 0 ) then
 			write(*,*) "### ERROR ### spinAvail.size() == 0", current, S, Si, Sj
 		end if
 
 		output = .true.
-
-		it1 => spinAvail.begin
-		do while( associated(it1) )
-
-			it2 => internalReactivesSpinAvail.begin
-			do while( associated(it2) )
-
-				if( abs( it1.data - it2.data ) < 0.1 ) then
+		do i=1,spinAvail.size()
+			do j=1,internalReactivesSpinAvail.size()
+				if( abs( spinAvail.at(i) - internalReactivesSpinAvail.at(j) ) < 0.1 ) then
 					output = .false.
 					return
 				end if
-
-				it2 => it2.next
 			end do
-
-			it1 => it1.next
 		end do
 
 		call spinAvail.clear()