hdiff output

r32049/sqnm.f90 2017-03-08 12:30:10.361011830 +0000 r32048/sqnm.f90 2017-03-08 12:30:10.621015301 +0000
347:                 atom1=BONDS(i,1)347:                 atom1=BONDS(i,1)
348:                 atom2=BONDS(i,2)348:                 atom2=BONDS(i,2)
349:                 !displacement is xyz(atom2)-xyz(atom1)349:                 !displacement is xyz(atom2)-xyz(atom1)
350:                 atom_disp(:)=xcoords(3*atom2-2:3*atom2)-xcoords(3*atom1-2:3*atom1)350:                 atom_disp(:)=xcoords(3*atom2-2:3*atom2)-xcoords(3*atom1-2:3*atom1)
351:                 bondvect(3*atom1-2:3*atom1,i)=atom_disp(:)351:                 bondvect(3*atom1-2:3*atom1,i)=atom_disp(:)
352:                 bondvect(3*atom2-2:3*atom2,i)=-atom_disp(:) !by antisymmetry352:                 bondvect(3*atom2-2:3*atom2,i)=-atom_disp(:) !by antisymmetry
353:             end do353:             end do
354:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then354:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then
355:                 print *, 'bond vectors'355:                 print *, 'bond vectors'
356:                 do i=1,min(writemax,nbonds) 356:                 do i=1,min(writemax,nbonds) 
357:                     write(*, '(f13.9)', advance='no') (bondvect(i,j),j=1,min(writemax,nbonds))357:                     write(*, '(*(f13.9))') (bondvect(i,j),j=1,min(writemax,nbonds))
358:                     write(*,*) "" 
359:                 end do358:                 end do
360:             end if 359:             end if 
361:             !FIND BOND VECTOR COEFFICIENTS 360:             !FIND BOND VECTOR COEFFICIENTS 
362:             !This is done by solving (bondvect^t * bondvect)*coeffs = bondvect^t * grad361:             !This is done by solving (bondvect^t * bondvect)*coeffs = bondvect^t * grad
363:             !find bondgrad=bondvect^t * grad and bondmat=(bondvect^t * bondvect)362:             !find bondgrad=bondvect^t * grad and bondmat=(bondvect^t * bondvect)
364:             do i=1,nbonds363:             do i=1,nbonds
365:                 bondcoeff(i)=dot(numcoords,bondvect(:,i),grad(:))364:                 bondcoeff(i)=dot(numcoords,bondvect(:,i),grad(:))
366:                 bondmat(i,i)=dot(numcoords,bondvect(:,i),bondvect(:,i))365:                 bondmat(i,i)=dot(numcoords,bondvect(:,i),bondvect(:,i))
367:                 do j=i+1,nbonds366:                 do j=i+1,nbonds
368:                     !compute only lower half since matrix is symmetric.367:                     !compute only lower half since matrix is symmetric.
374:             do i=1,nbonds373:             do i=1,nbonds
375:                 if ((oldcoeff(i)>=0.0d0 .and. bondcoeff(i)>=0.0d0) .or. (oldcoeff(i)<0.0d0 .and. bondcoeff(i)<0.0d0)) then374:                 if ((oldcoeff(i)>=0.0d0 .and. bondcoeff(i)>=0.0d0) .or. (oldcoeff(i)<0.0d0 .and. bondcoeff(i)<0.0d0)) then
376:                     signflip=signflip+1375:                     signflip=signflip+1
377:                 end if 376:                 end if 
378:             end do377:             end do
379:             oldcoeff(:)=bondcoeff(:)378:             oldcoeff(:)=bondcoeff(:)
380: 379: 
381:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then380:             if(debug_run .and. SQNM_DEBUGLEVEL>=1) then
382:                 print *, 'bond matrix'381:                 print *, 'bond matrix'
383:                 do i=1,min(writemax,nbonds) 382:                 do i=1,min(writemax,nbonds) 
384:                     write(*, '(f13.9)', advance='no') (bondmat(i,j),j=1,min(writemax,nbonds))383:                     write(*, '(*(f13.9))') (bondmat(i,j),j=1,min(writemax,nbonds))
385:                     write(*,*) "" 
386:                 end do384:                 end do
387:                 print *, 'bond coeffs'385:                 print *, 'bond coeffs'
388:                 do i=1,min(writemax,nbonds)386:                 do i=1,min(writemax,nbonds)
389:                     write(*,'(f16.10)') bondcoeff(i)387:                     write(*,'(f16.10)') bondcoeff(i)
390:                 end do388:                 end do
391:                 print *, 'Number of sign flips: ', signflip, '; alpha_bio = ', alpha_bio389:                 print *, 'Number of sign flips: ', signflip, '; alpha_bio = ', alpha_bio
392:             end if 390:             end if 
393: 391: 
394: 392: 
395:             !Solve linear system Bondmat*coeff=bondvect. 393:             !Solve linear system Bondmat*coeff=bondvect. 
446:             graddisp(:,1)=grad(:)-gradhist(:,1)444:             graddisp(:,1)=grad(:)-gradhist(:,1)
447:             !for remaining columns, just use historical positions.445:             !for remaining columns, just use historical positions.
448:             do i=2,nhist446:             do i=2,nhist
449:                 xdisp(:,i)=xhist(:,i-1)-xhist(:,i)447:                 xdisp(:,i)=xhist(:,i-1)-xhist(:,i)
450:                 graddisp(:,i)=gradhist(:,i-1)-gradhist(:,i)448:                 graddisp(:,i)=gradhist(:,i-1)-gradhist(:,i)
451:             end do  449:             end do  
452: 450: 
453:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then 451:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then 
454:                 write(*, '(A)') 'xhist'452:                 write(*, '(A)') 'xhist'
455:                 do i=1,min(writemax,numcoords)453:                 do i=1,min(writemax,numcoords)
456:                     write(*, '(f20.15)', advance='no') (xhist(i,j),j=1,nhist)454:                     write(*, '(*(f20.15))') (xhist(i,j),j=1,nhist)
457:                     write(*,*) "" 
458:                 end do 455:                 end do 
459:                 write(*, '(A)') 'gradhist'456:                 write(*, '(A)') 'gradhist'
460:                 do i=1,min(writemax,numcoords)457:                 do i=1,min(writemax,numcoords)
461:                     write(*, '(f20.15)', advance='no') (gradhist(i,j),j=1,nhist)458:                     write(*, '(*(f20.15))') (gradhist(i,j),j=1,nhist)
462:                     write(*,*) "" 
463:                 end do 459:                 end do 
464:                 write(*, '(A)') 'xdisp'460:                 write(*, '(A)') 'xdisp'
465:                 do i=1,min(writemax,numcoords)461:                 do i=1,min(writemax,numcoords)
466:                     write(*, '(f16.10)', advance='no') (xdisp(i,j),j=1,nhist)462:                     write(*, '(*(f16.10))') (xdisp(i,j),j=1,nhist)
467:                     write(*,*) "" 
468:                 end do463:                 end do
469:                 write(*, '(A)') 'graddisp'464:                 write(*, '(A)') 'graddisp'
470:                 do i=1,min(writemax,numcoords)465:                 do i=1,min(writemax,numcoords)
471:                     write(*, '(f16.10)', advance='no') (graddisp(i,j),j=1,nhist)466:                     write(*, '(*(f16.10))') (graddisp(i,j),j=1,nhist)
472:                     write(*,*) "" 
473:                 end do   467:                 end do   
474:             end if 468:             end if 
475: 469: 
476:             !FIND NORMS AND NORMALIZE X-DISPLACEMENT VECTORS470:             !FIND NORMS AND NORMALIZE X-DISPLACEMENT VECTORS
477:             allocate(norms(nhist))471:             allocate(norms(nhist))
478:             do i=1, nhist 472:             do i=1, nhist 
479:                 norms(i)=norm(numcoords,xdisp(:,i))473:                 norms(i)=norm(numcoords,xdisp(:,i))
480:                 if (norms(i)<1.0d-100) then !floating point errors caused underflow. Return to avoid division by 0.474:                 if (norms(i)<1.0d-100) then !floating point errors caused underflow. Return to avoid division by 0.
481:                     if (debug_run  .and. SQNM_DEBUGLEVEL>=0) print *,'Norm too small for division on run ',runs,', iter,',iter,&475:                     if (debug_run  .and. SQNM_DEBUGLEVEL>=0) print *,'Norm too small for division on run ',runs,', iter,',iter,&
482:                         'hist ', nhist, ' Bombing out...'476:                         'hist ', nhist, ' Bombing out...'
483:                     converget=.false.477:                     converget=.false.
484:                     return478:                     return
485:                 end if479:                 end if
486:                 !non-normalized x-displacements will not be needed again, so we can normalize in original matrix.480:                 !non-normalized x-displacements will not be needed again, so we can normalize in original matrix.
487:                 xdisp(:,i)=xdisp(:,i)/norms(i)481:                 xdisp(:,i)=xdisp(:,i)/norms(i)
488:             end do482:             end do
489: 483: 
490:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then 484:             if (debug_run .and. SQNM_DEBUGLEVEL>=2) then 
491:                 write(*, '(A)') 'displacement norms'485:                 write(*, '(A)') 'displacement norms'
492:                 do i=1,min(writemax,nhist)486:                 do i=1,min(writemax,nhist)
493:                     write(*, '(f16.10)', advance='no') (norms(i))487:                     write(*, '(*(f16.10))') (norms(i))
494:                     write(*,*) "" 
495:                 end do 488:                 end do 
496:                 write(*, '(A)') 'xdisp - normalized'489:                 write(*, '(A)') 'xdisp - normalized'
497:                 do i=1,min(writemax,numcoords)490:                 do i=1,min(writemax,numcoords)
498:                     write(*, '(f16.10)', advance='no') (xdisp(i,j),j=1,min(writemax,nhist))491:                     write(*, '(*(f16.10))') (xdisp(i,j),j=1,min(writemax,nhist))
499:                     write(*,*) "" 
500:                 end do  492:                 end do  
501:             end if493:             end if
502: 494: 
503:             !FIND OVERLAP MATRIX495:             !FIND OVERLAP MATRIX
504:             allocate(overlap(nhist, nhist))496:             allocate(overlap(nhist, nhist))
505:             do i=1, nhist497:             do i=1, nhist
506:                 !diagonal elements498:                 !diagonal elements
507:                 overlap(i,i)=dot(numcoords,xdisp(:,i),xdisp(:,i))499:                 overlap(i,i)=dot(numcoords,xdisp(:,i),xdisp(:,i))
508:                 do j=i+1, nhist500:                 do j=i+1, nhist
509:                     !compute only lower half since matrix is symmetric.501:                     !compute only lower half since matrix is symmetric.
510:                     overlap(j,i)=dot(numcoords,xdisp(:,i),xdisp(:,j))502:                     overlap(j,i)=dot(numcoords,xdisp(:,i),xdisp(:,j))
511:                 end do503:                 end do
512:             end do504:             end do
513: 505: 
514:             if (debug_run .and. SQNM_DEBUGLEVEL>=3) then 506:             if (debug_run .and. SQNM_DEBUGLEVEL>=3) then 
515:                 write(*, '(A)') 'overlap'507:                 write(*, '(A)') 'overlap'
516:                 do i=1,nhist508:                 do i=1,nhist
517:                     write(*, '(f16.10)', advance='no') (overlap(i,j),j=1,nhist)509:                     write(*, '(*(f16.10))') (overlap(i,j),j=1,nhist)
518:                     write(*,*) "" 
519:                 end do510:                 end do
520:             end if511:             end if
521: 512: 
522:             !GET EIGENVALUES AND EIGENVECTORS OF OVERLAP MATRIX513:             !GET EIGENVALUES AND EIGENVECTORS OF OVERLAP MATRIX
523:             !allocate vectors for eigenvectors and for work for dsyev514:             !allocate vectors for eigenvectors and for work for dsyev
524:             allocate(eigenvalues(nhist))515:             allocate(eigenvalues(nhist))
525:             !dsyev will place eigenvalues of overlap into the eigenvalues vector in order of ascending value, 516:             !dsyev will place eigenvalues of overlap into the eigenvalues vector in order of ascending value, 
526:             !and the columns of overlap will be replaces with the corresponding eigenvectors.517:             !and the columns of overlap will be replaces with the corresponding eigenvectors.
527:             !Note that the LAPACK function for symmetric matrices is called since we calculated the lower518:             !Note that the LAPACK function for symmetric matrices is called since we calculated the lower
528:             !triangular portion of overlap.519:             !triangular portion of overlap.
529:             call dsyev('V',"L",nhist,overlap,nhist,eigenvalues,work,lwork,infoc) 520:             call dsyev('V',"L",nhist,overlap,nhist,eigenvalues,work,lwork,infoc) 
530:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 521:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 
531:                 write(*, '(A)') 'eigenvalues'522:                 write(*, '(A)') 'eigenvalues'
532:                 do i=1,nhist523:                 do i=1,nhist
533:                     write(*, '(f16.10)', advance='no') (eigenvalues(i))524:                     write(*, '(*(f16.10))') (eigenvalues(i))
534:                     write(*,*) "" 
535:                 end do525:                 end do
536:                 write(*, '(A)') 'overlap matrix eigenvectors'526:                 write(*, '(A)') 'overlap matrix eigenvectors'
537:                 do i=1,nhist527:                 do i=1,nhist
538:                     write(*, '(f16.10)', advance='no') (overlap(i,j),j=1,nhist)528:                     write(*, '(*(f16.10))') (overlap(i,j),j=1,nhist)
539:                     write(*,*) ""529:                 end do   
540:                 end do   
541:             end if 530:             end if 
542: 531: 
543:             !FIND THE DIMENSIONS OF THE SIGNIFICANT SUBSPACE. This is done by comparing the ratio of all eigenvalues to 532:             !FIND THE DIMENSIONS OF THE SIGNIFICANT SUBSPACE. This is done by comparing the ratio of all eigenvalues to 
544:             !the largest eigenvalue. This is used to find the dimension of the significant subspace:533:             !the largest eigenvalue. This is used to find the dimension of the significant subspace:
545:             ndim=1 !minimum dimension if nhist>0534:             ndim=1 !minimum dimension if nhist>0
546:             maxlambda=1.0D0/eigenvalues(nhist) !eigenvalues are in ascending order, so eigenvalues(nhist) is the largest.535:             maxlambda=1.0D0/eigenvalues(nhist) !eigenvalues are in ascending order, so eigenvalues(nhist) is the largest.
547:             do i=1, nhist-1536:             do i=1, nhist-1
548:                 if (maxlambda*eigenvalues(i)>=cutoffratio) then !equivalent to checking eigenvalues(i)/eigenvalues(nhist)537:                 if (maxlambda*eigenvalues(i)>=cutoffratio) then !equivalent to checking eigenvalues(i)/eigenvalues(nhist)
549:                     ndim=ndim+1538:                     ndim=ndim+1
550:                 end if 539:                 end if 
563:                 !note that any eigenvalues should be at least zero since the overlap matrix is positive semi-definite; 552:                 !note that any eigenvalues should be at least zero since the overlap matrix is positive semi-definite; 
564:                 !however, due to floating point error, we may have very slightly negative or zero eigenvalues. Take max(). 553:                 !however, due to floating point error, we may have very slightly negative or zero eigenvalues. Take max(). 
565:                 tempval=1.0D0/sqrt(max(eigenvalues(i),1.0d-10))554:                 tempval=1.0D0/sqrt(max(eigenvalues(i),1.0d-10))
566:                 subspan(:,i)=subspan(:,i)*tempval555:                 subspan(:,i)=subspan(:,i)*tempval
567:                 gradspan(:,i)=gradspan(:,i)*tempval556:                 gradspan(:,i)=gradspan(:,i)*tempval
568:             end do557:             end do
569: 558: 
570:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 559:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 
571:                 write(*, '(A)') 'subspance span'560:                 write(*, '(A)') 'subspance span'
572:                 do i=1,min(writemax,numcoords)561:                 do i=1,min(writemax,numcoords)
573:                    write(*, '(f16.10)', advance='no') (subspan(i,j),j=1,min(writemax,nhist))562:                     write(*, '(*(f16.10))') (subspan(i,j),j=1,min(writemax,nhist))
574:                     write(*,*) "" 
575:                 end do 563:                 end do 
576:                 write(*, '(A)') 'subspance gradients'564:                 write(*, '(A)') 'subspance gradients'
577:                 do i=1,min(writemax,numcoords)565:                 do i=1,min(writemax,numcoords)
578:                     write(*, '(f16.10)', advance='no') (gradspan(i,j),j=1,min(writemax,nhist))566:                     write(*, '(*(f16.10))') (gradspan(i,j),j=1,min(writemax,nhist))
579:                     write(*,*) "" 
580:                 end do 567:                 end do 
581:             end if 568:             end if 
582: 569: 
583:             !SELECT VECTORS AND GRADIENTS IN THE SIGNIFICANT SUBSPACE570:             !SELECT VECTORS AND GRADIENTS IN THE SIGNIFICANT SUBSPACE
584:             !there will be ndim of these vectors, out of nhist, with ndim<=nhist571:             !there will be ndim of these vectors, out of nhist, with ndim<=nhist
585:             allocate(sigsubspan(numcoords,ndim))572:             allocate(sigsubspan(numcoords,ndim))
586:             allocate(siggradspan(numcoords,ndim))573:             allocate(siggradspan(numcoords,ndim))
587:             !Note that eigenvalues were listed in ascending order before, so we reverse the direction of vector selection.574:             !Note that eigenvalues were listed in ascending order before, so we reverse the direction of vector selection.
588:             !The first vectors correspond with the largest eigenvalue and so on.575:             !The first vectors correspond with the largest eigenvalue and so on.
589:             do i=1, ndim576:             do i=1, ndim
590:                 sigsubspan(:,i)=subspan(:,nhist-i+1)577:                 sigsubspan(:,i)=subspan(:,nhist-i+1)
591:                 siggradspan(:,i)=gradspan(:,nhist-i+1)578:                 siggradspan(:,i)=gradspan(:,nhist-i+1)
592:             end do579:             end do
593: 580: 
594:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 581:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 
595:                 write(*, '(A)') 'significant subspance span'582:                 write(*, '(A)') 'significant subspance span'
596:                 do i=1,min(writemax,numcoords)583:                 do i=1,min(writemax,numcoords)
597:                     write(*, '(f16.10)', advance='no') (sigsubspan(i,j),j=1,min(writemax,ndim))584:                     write(*, '(*(f16.10))') (sigsubspan(i,j),j=1,min(writemax,ndim))
598:                     write(*,*) "" 
599:                 end do 585:                 end do 
600:                 write(*, '(A)') 'significant subspance gradients'586:                 write(*, '(A)') 'significant subspance gradients'
601:                 do i=1,min(writemax,numcoords)587:                 do i=1,min(writemax,numcoords)
602:                     write(*, '(f16.10)', advance='no') (siggradspan(i,j),j=1,min(writemax,ndim))588:                     write(*, '(*(f16.10))') (siggradspan(i,j),j=1,min(writemax,ndim))
603:                     write(*,*) "" 
604:                 end do 589:                 end do 
605:             end if 590:             end if 
606: 591: 
607:             !PRODUCE SYMMETRIC SHAPE OPERATOR592:             !PRODUCE SYMMETRIC SHAPE OPERATOR
608:             allocate(SO(ndim,ndim))593:             allocate(SO(ndim,ndim))
609:             SO(:,:)=0.0D0594:             SO(:,:)=0.0D0
610:             do i=1, ndim595:             do i=1, ndim
611:                 SO(i,i)=dot(numcoords,siggradspan(:,i),sigsubspan(:,i)) !diagonals of shape operator596:                 SO(i,i)=dot(numcoords,siggradspan(:,i),sigsubspan(:,i)) !diagonals of shape operator
612:                 do j=i+1, ndim597:                 do j=i+1, ndim
613:                     !matrix is symmetric; find only lower triangular 598:                     !matrix is symmetric; find only lower triangular 
614:                     SO(j,i)=0.5D0*(dot(numcoords,siggradspan(:,i),sigsubspan(:,j))+ &599:                     SO(j,i)=0.5D0*(dot(numcoords,siggradspan(:,i),sigsubspan(:,j))+ &
615:                         dot(numcoords,siggradspan(:,j),sigsubspan(:,i))) 600:                         dot(numcoords,siggradspan(:,j),sigsubspan(:,i))) 
616:                 end do601:                 end do
617:             end do602:             end do
618: 603: 
619:             if (debug_run .and. SQNM_DEBUGLEVEL>=5) then 604:             if (debug_run .and. SQNM_DEBUGLEVEL>=5) then 
620:                 write(*, '(A)') 'shape operator matrix'605:                 write(*, '(A)') 'shape operator matrix'
621:                 do i=1,min(writemax,ndim)606:                 do i=1,min(writemax,ndim)
622:                     write(*, '(f16.10)', advance='no') (SO(i,j),j=1,min(writemax,ndim))607:                     write(*, '(*(f16.10))') (SO(i,j),j=1,min(writemax,ndim))
623:                     write(*,*) "" 
624:                 end do   608:                 end do   
625:             end if 609:             end if 
626: 610: 
627:             !PERFORM EIGENVALUE DECOMFC=ifort cmake -DCMAKE_BUILD_TYPE=Debug -DWITH_AMBER=1 ../../sourcePOSITION OF SHAPE OPERATOR611:             !PERFORM EIGENVALUE DECOMFC=ifort cmake -DCMAKE_BUILD_TYPE=Debug -DWITH_AMBER=1 ../../sourcePOSITION OF SHAPE OPERATOR
628:             !we no longer need the eigenvectors from the overlap matrix, so we reuse the vector.612:             !we no longer need the eigenvectors from the overlap matrix, so we reuse the vector.
629:             deallocate(eigenvalues)613:             deallocate(eigenvalues)
630:             allocate(eigenvalues(ndim))614:             allocate(eigenvalues(ndim))
631:             !Note that the LAPACK function for symmetric matrices is called since we calculated the lower615:             !Note that the LAPACK function for symmetric matrices is called since we calculated the lower
632:             !triangular portion of the shape operator.616:             !triangular portion of the shape operator.
633:             call dsyev('V',"L",ndim,SO,ndim,eigenvalues,work,lwork,infoc)  617:             call dsyev('V',"L",ndim,SO,ndim,eigenvalues,work,lwork,infoc)  
634: 618: 
635:             if (debug_run .and. SQNM_DEBUGLEVEL>=5) then 619:             if (debug_run .and. SQNM_DEBUGLEVEL>=5) then 
636:                 write(*, '(A)') 'shape operator eigenvalues'620:                 write(*, '(A)') 'shape operator eigenvalues'
637:                 do i=1,min(writemax,ndim)621:                 do i=1,min(writemax,ndim)
638:                     write(*, '(f16.10)') (eigenvalues(i))622:                     write(*, '(f16.10)') (eigenvalues(i))
639:                 end do623:                 end do
640:                 write(*, '(A)') 'shape operator matrix eigenvectors'624:                 write(*, '(A)') 'shape operator matrix eigenvectors'
641:                 do i=1,min(writemax,ndim)625:                 do i=1,min(writemax,ndim)
642:                     write(*, '(f16.10)', advance='no') (SO(i,j),j=1,min(writemax,ndim))626:                     write(*, '(*(f16.10))') (SO(i,j),j=1,min(writemax,ndim))
643:                     write(*,*) "" 
644:                 end do   627:                 end do   
645:             end if 628:             end if 
646: 629: 
647:             !FIND PROJECTED PRINCIPAL DIRECTIONS630:             !FIND PROJECTED PRINCIPAL DIRECTIONS
648:             allocate(principals(numcoords,ndim))631:             allocate(principals(numcoords,ndim))
649:             do i=1, ndim632:             do i=1, ndim
650:                 principals(:,i)=0.0D0633:                 principals(:,i)=0.0D0
651:                 do j=1, ndim634:                 do j=1, ndim
652:                     principals(:,i)=principals(:,i)+SO(j,i)*sigsubspan(:,j)635:                     principals(:,i)=principals(:,i)+SO(j,i)*sigsubspan(:,j)
653:                 end do636:                 end do
654:             end do 637:             end do 
655: 638: 
656:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 639:             if (debug_run .and. SQNM_DEBUGLEVEL>=4) then 
657:                 write(*, '(A)') 'principals'640:                 write(*, '(A)') 'principals'
658:                 do i=1,min(writemax,numcoords)641:                 do i=1,min(writemax,numcoords)
659:                     write(*, '(f16.10)', advance='no') (principals(i,j),j=1,min(writemax,ndim))642:                     write(*, '(*(f16.10))') (principals(i,j),j=1,min(writemax,ndim))
660:                     write(*,*) "" 
661:                 end do  643:                 end do  
662:             end if644:             end if
663: 645: 
664:             !ADJUST CURVATURES BY FINDING RESIDUES646:             !ADJUST CURVATURES BY FINDING RESIDUES
665:             allocate(residues(ndim))647:             allocate(residues(ndim))
666:             !calculate residues648:             !calculate residues
667:             do i=1, ndim649:             do i=1, ndim
668:                 residues(i)=0.0D0650:                 residues(i)=0.0D0
669:                 residuetemp(:)=0.0D0651:                 residuetemp(:)=0.0D0
670:                 do j =1, ndim652:                 do j =1, ndim
688: 670: 
689:             !FORM SUBSPACE PORTION OF PRECONDITIONED GRADIENT671:             !FORM SUBSPACE PORTION OF PRECONDITIONED GRADIENT
690:             pregrad(:)=0.0D0672:             pregrad(:)=0.0D0
691:             do i=1, ndim673:             do i=1, ndim
692:                 tempval=dot(numcoords,grad,principals(:,i))/eigenvalues(i)674:                 tempval=dot(numcoords,grad,principals(:,i))/eigenvalues(i)
693:                 pregrad(:)=pregrad(:)+tempval*principals(:,i)675:                 pregrad(:)=pregrad(:)+tempval*principals(:,i)
694:             end do676:             end do
695:             if (debug_run .and. SQNM_DEBUGLEVEL>=3) then 677:             if (debug_run .and. SQNM_DEBUGLEVEL>=3) then 
696:                 write(*, '(A)') 'Preconditioned Gradient in Subspace'678:                 write(*, '(A)') 'Preconditioned Gradient in Subspace'
697:                 do i=1,min(writemax,numcoords)679:                 do i=1,min(writemax,numcoords)
698:                     write(*, '(f16.10)', advance='no') pregrad(i)680:                     write(*, '(*(f16.10))') pregrad(i)
699:                     write(*,*) "" 
700:                 end do681:                 end do
701:             end if682:             end if
702:             !FIND COS(ANGLE) BETWEEN GRADIENT AND PRECONDITIONED GRADIENT TO ADJUST ALPHA LATER.683:             !FIND COS(ANGLE) BETWEEN GRADIENT AND PRECONDITIONED GRADIENT TO ADJUST ALPHA LATER.
703:             !ALSO FIND ORTHOGONAL COMPLEMENT OF GRADIENT: gradperp=(I-P')grad; invproj=(I-P')684:             !ALSO FIND ORTHOGONAL COMPLEMENT OF GRADIENT: gradperp=(I-P')grad; invproj=(I-P')
704:             angle=dot(numcoords,pregrad,grad)/(norm(numcoords,pregrad)*norm(numcoords,grad))685:             angle=dot(numcoords,pregrad,grad)/(norm(numcoords,pregrad)*norm(numcoords,grad))
705:             !Note P'=principals*principals^T; can use dot products for each entry of invproj; add 1 to diagonals to get (I-P')686:             !Note P'=principals*principals^T; can use dot products for each entry of invproj; add 1 to diagonals to get (I-P')
706:             do i=1, numcoords687:             do i=1, numcoords
707:                 invproj(i,i)=1.0d0-dot(ndim,principals(i,:),principals(i,:))688:                 invproj(i,i)=1.0d0-dot(ndim,principals(i,:),principals(i,:))
708:                 do j=i+1, numcoords689:                 do j=i+1, numcoords
709:                     invproj(j,i)=-dot(ndim,principals(i,:),principals(j,:))690:                     invproj(j,i)=-dot(ndim,principals(i,:),principals(j,:))


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0