assign_strbnd Subroutine

    subroutine assign_strbnd(bds, prm_buf)
        use mod_memory, only: mallocate, mfree
        use mod_bonded, only: strbnd_init 
        use mod_constants, only: kcalmol2au, angstrom2au
        
        implicit none
       
        type(ommp_bonded_type), intent(inout) :: bds
        !! Bonded potential data structure
        character(len=OMMP_STR_CHAR_MAX), intent(in) :: prm_buf(:)
        !! Char buffer containing the prm file loaded in RAM

        integer(ip) :: il, i, j, tokb, toke, isb, nstrbnd, &
                       cla, clb, clc, a, b, c, jc, jb, maxsb, &
                       l1a, l1b, l2a, l2b
        character(len=OMMP_STR_CHAR_MAX) :: line, errstring
        integer(ip), allocatable :: classa(:), classb(:), classc(:), sbtmp(:), &
                                    sbattmp(:, :), at2bnd(:), at2ang(:)
        real(rp), allocatable :: k1(:), k2(:)
        logical :: done, thet_done, l1_done, l2_done
        type(ommp_topology_type), pointer :: top

        top => bds%top

        if(.not. top%atclass_initialized .or. .not. top%atz_initialized) then
            call read_atom_cards(top, prm_buf)
        end if
        
        ! Read all the lines of file just to count how large vector should be 
        ! allocated
        nstrbnd = 0
        do il=1, size(prm_buf) 
            line = prm_buf(il)
            if(line(:7) == 'strbnd ') nstrbnd = nstrbnd + 1
        end do

        maxsb = (top%conn(2)%ri(top%mm_atoms+1)-1) / 2
        call mallocate('assign_strbnd [classa]', nstrbnd, classa)
        call mallocate('assign_strbnd [classb]', nstrbnd, classb)
        call mallocate('assign_strbnd [classc]', nstrbnd, classc)
        call mallocate('assign_strbnd [eqang]', nstrbnd, k1)
        call mallocate('assign_strbnd [kang]', nstrbnd, k2)
        call mallocate('assign_strbnd [sbtmp]', maxsb, sbtmp)
        call mallocate('assign_strbnd [sbattmp]', 3_ip, maxsb, sbattmp)
        call mallocate('assign_strbnd [at2bnd]', top%mm_atoms, at2bnd)
        call mallocate('assign_strbnd [at2ang]', top%mm_atoms, at2ang)

        ! Restart the reading from the beginning to actually save the parameters
        isb = 1
        i=1
        do il=1, size(prm_buf) 
            line = prm_buf(il)
           
            if(line(:7) == 'strbnd ') then
                tokb = 8
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong STRBND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) classa(isb)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong STRBND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) classb(isb)
                
                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isint(line(tokb:toke))) then
                    write(errstring, *) "Wrong STRBND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) classc(isb)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong STRBND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) k1(isb)

                tokb = toke + 1
                toke = tokenize(line, tokb)
                if(.not. isreal(line(tokb:toke))) then
                    write(errstring, *) "Wrong STRBND card"
                    call fatal_error(errstring)
                end if
                read(line(tokb:toke), *) k2(isb)

                isb = isb + 1
            end if
            i = i+1
        end do
        
        isb = 1
        do a=1, top%mm_atoms
            cla = top%atclass(a)
            do jb=top%conn(2)%ri(a), top%conn(2)%ri(a+1)-1
                b = top%conn(2)%ci(jb)
                if(a > b) cycle
                clb = top%atclass(b)
                
                do jc=top%conn(1)%ri(a), top%conn(1)%ri(a+1)-1
                    c = top%conn(1)%ci(jc)
                    if(all(top%conn(1)%ci(top%conn(1)%ri(b):&
                                          top%conn(1)%ri(b+1)-1) /= c)) cycle
                    ! There is an angle in the form A-C-B
                    clc = top%atclass(c)
                    done = .false.

                    do j=1, nstrbnd
                        if((cla == classa(j) &
                            .and. clb == classc(j) &
                            .and. clc == classb(j)) .or. &
                           (clb == classa(j) &
                            .and. cla == classc(j) &
                            .and. clc == classb(j))) then
                            sbattmp(1,isb) = a
                            sbattmp(2,isb) = c
                            sbattmp(3,isb) = b
                            if(cla == classa(j)) then
                                ! Assign the correct k to each bond stretching!
                                sbtmp(isb) = j
                            else
                                sbtmp(isb) = -j
                            end if
                            isb = isb + 1
                            exit
                        end if
                    end do
                end do
            end do
        end do

        call strbnd_init(bds, isb-1)
        if(isb-1 < 1) then
            !! No parameters are defined, nothing to do.
            return
        end if

        at2bnd(1) = 1
        do i=1, top%mm_atoms-1
            at2bnd(i+1) = at2bnd(i)
            do j=top%conn(1)%ri(i), top%conn(1)%ri(i+1)-1
                if(i < top%conn(1)%ci(j)) then
                    at2bnd(i+1) = at2bnd(i+1) + 1
                end if
            end do
        end do
            
        at2ang = 0
        do j=1, size(bds%angleat, 2)
          if(at2ang(bds%angleat(1,j)) == 0) at2ang(bds%angleat(1,j)) = j
        end do

        !$omp parallel do default(shared) schedule(dynamic) &
        !$omp private(i,j,l1a,l1b,l2a,l2b,l1_done,l2_done,thet_done)
        do i=1, isb-1
            ! First assign the parameters
            bds%strbndat(:,i) = sbattmp(:,i)
            j = abs(sbtmp(i))
            if(sbtmp(i) > 0) then
                bds%strbndk1(i) = k1(j) * kcalmol2au / angstrom2au
                bds%strbndk2(i) = k2(j) * kcalmol2au / angstrom2au
            else
                bds%strbndk1(i) = k2(j) * kcalmol2au / angstrom2au
                bds%strbndk2(i) = k1(j) * kcalmol2au / angstrom2au
            end if
            
            l1a = min(bds%strbndat(1,i), bds%strbndat(2,i))
            l1b = max(bds%strbndat(1,i), bds%strbndat(2,i))
            l2a = min(bds%strbndat(3,i), bds%strbndat(2,i))
            l2b = max(bds%strbndat(3,i), bds%strbndat(2,i))
           
            ! Now search for the corresponding bond and angle parameters to
            ! set the equilibrium distances and angle
            l1_done = .false.
            l2_done = .false.
            thet_done = .false.

            do j=at2bnd(l1a), size(bds%bondat, 2)
                if(l1a == bds%bondat(1,j) .and. l1b == bds%bondat(2,j)) then
                    l1_done = .true.
                    bds%strbndl10(i) = bds%l0bond(j)
                    exit
                end if
            end do

            do j=at2bnd(l2a), size(bds%bondat, 2)
                if(l2a == bds%bondat(1,j) .and. l2b == bds%bondat(2,j)) then
                    l2_done = .true.
                    bds%strbndl20(i) = bds%l0bond(j)
                    exit
                end if
            end do
            
            do j=at2ang(bds%strbndat(1,i)), size(bds%angleat, 2)
                if(all(bds%strbndat(:,i) == bds%angleat(:,j))) then
                    thet_done = .true.
                    bds%strbndthet0(i) = bds%eqangle(j)
                    exit
                end if
            end do
            
            if(.not.(l1_done .and. l2_done .and. thet_done)) then
                call fatal_error("Ill-defined stretch-bending cross term")
            end if
        end do

        call mfree('assign_strbnd [classa]', classa)
        call mfree('assign_strbnd [classb]', classb)
        call mfree('assign_strbnd [classc]', classc)
        call mfree('assign_strbnd [eqang]', k1)
        call mfree('assign_strbnd [kang]', k2)
        call mfree('assign_strbnd [sbtmp]', sbtmp)
        call mfree('assign_strbnd [sbattmp]', sbattmp)
        call mfree('assign_strbnd [at2bnd]', at2bnd)
        call mfree('assign_strbnd [at2ang]', at2ang)
    
    end subroutine assign_strbnd