-
Notifications
You must be signed in to change notification settings - Fork 149
/
gfnff_ini.f90
2032 lines (1896 loc) · 89.8 KB
/
gfnff_ini.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! This file is part of xtb.
!
! Copyright (C) 2019-2020 Stefan Grimme
!
! xtb is free software: you can redistribute it and/or modify it under
! the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! xtb is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see <https://www.gnu.org/licenses/>.
module xtb_gfnff_ini
contains
subroutine gfnff_ini(env,pr,makeneighbor,mol,gen,param,topo,accuracy)
use xtb_mctc_accuracy, only : wp, sp
use xtb_type_molecule
use xtb_type_environment, only : TEnvironment
use xtb_gfnff_param, only : efield, gfnff_thresholds
use xtb_gfnff_data, only : TGFFData
use xtb_gfnff_topology, only : TGFFTopology
use xtb_gfnff_generator, only : TGFFGenerator
use xtb_gfnff_ini2
use xtb_gfnff_eg, only : gfnff_dlogcoord
use xtb_gfnff_mrec, only : mrecgff
use xtb_gfnff_fraghess
use xtb_restart
use xtb_mctc_constants
implicit none
character(len=*), parameter :: source = 'gfnff_ini'
!--------------------------------------------------------------------------------------------------
type(TEnvironment), intent(inout) :: env
type(TMolecule), intent(in) :: mol ! # molecule type
type(TGFFTopology), intent(inout) :: topo
type(TGFFGenerator), intent(in) :: gen
type(TGFFData), intent(in) :: param
real(wp), intent(in) :: accuracy
logical, intent(in) :: pr ! print flag
logical, intent(in) :: makeneighbor ! make a neigbor list or use existing one?
!--------------------------------------------------------------------------------------------------
integer ati,atj,atk,i,j,k,l,lin,nn,ii,jj,kk,ll,m,rings,ia,ja,ij,ix,nnn,idum,ip,ji,no,nbi
integer ineig,jneig,nrot,bbtyp,ringtyp,nn1,nn2,hybi,hybj,pis,ka,nh,jdum,hcalc,nc
integer ringsi,ringsj,ringsk,ringl,npi,nelpi,picount,npiall,maxtors,rings4,nheav
integer nm,maxhb,ki,n13,current,ncarbo,mtyp1,mtyp2
integer ind3(3),sr(20),cr(10,20),niel(86)
integer qloop_count,nf,nsi,nmet,nhi,nhj,ifrag
integer hbA,hbH,Bat,atB,Aat,Hat
integer AHB_nr
integer bond_hbn
interface
integer function itabrow6(i)
integer i
end function
end interface
real(wp) r0,ff,omega,f1,f2,phi,valijklff,ringf,fcn
real(wp) shift,dum,dum1,dum2,dum4,qafac,fqq,feta
real(wp) sumppi,fpi,fxh,fijk,fsrb2,ees
real(wp) fheavy,fn,eold,fctot,fij
real(wp) hbpi(2),hbpj(2),sdum3(3)
real(wp) bstrength
real(wp) xx(20)
real(wp) fkl,qreps,fbsmall,bohr
parameter (bohr=1.0_wp/0.52917726_wp)
real(wp), parameter :: rabd_cutoff = 13.0_wp
logical lring,picon,notpicon,bridge,sp3ij,ccij,success
logical heavy,triple,piat,sp3kl,ex,cnij,frag_charges_known
integer,allocatable :: btyp(:),imetal(:),nbm(:,:),nbf(:,:)
integer,allocatable :: itag(:)
integer,allocatable :: piadr(:),piadr2(:),piadr3(:),piadr4(:)
integer,allocatable :: itmp(:),sring(:,:),cring(:,:,:)
integer,allocatable :: ipis(:),pimvec(:),nbpi(:,:),piel(:)
integer,allocatable :: lin_AHB(:)
integer,allocatable :: bond_hbl(:,:)
real(wp),allocatable:: rab (:)
real(wp),allocatable:: sqrab(:)
real(wp),allocatable:: cn (:)
real(wp),allocatable:: dcn(:,:,:)
real(wp),allocatable:: dgam(:), dxi(:)
real(wp),allocatable:: mchar(:)
real(wp),allocatable:: rtmp (:)
real(wp),allocatable:: pbo (:)
real(wp),allocatable:: qtmp (:), dqa(:), qah(:)
real(wp),allocatable:: Api(:,:),S(:,:),Pold(:,:),pibo(:),occ(:),eps(:)
real(wp),allocatable:: pispop(:),pisea(:),pisip(:),apisave(:,:)
real(sp),allocatable:: rabd(:,:)
character(len=255) atmp
integer :: ich, err
real(wp) :: dispthr, cnthr, repthr, hbthr1, hbthr2
logical :: exitRun
call gfnff_thresholds(accuracy, dispthr, cnthr, repthr, hbthr1, hbthr2)
if (pr) then
write(env%unit,*)
write(env%unit,'(10x,"entering GFN-FF setup routine... ",i0)') mol%n
endif
write(env%unit,*)
write(env%unit,'(10x,"==================== Thresholds ====================")')
write(env%unit,'(10x,"CN :",f12.5)') cnthr
write(env%unit,'(10x,"rep :",f12.5)') repthr
write(env%unit,'(10x,"disp:",f12.5)') dispthr
write(env%unit,'(10x,"HB1 :",f12.5)') hbthr1
write(env%unit,'(10x,"HB2 :",f12.5)') hbthr2
write(env%unit,*)
allocate( rab(mol%n*(mol%n+1)/2), source = 0.0d0 )
allocate( cn(mol%n), source = 0.0d0 )
allocate( sqrab(mol%n*(mol%n+1)/2), source = 0.0d0 )
allocate( topo%hyb(mol%n), source = 0 )
allocate( topo%alphanb(mol%n*(mol%n+1)/2), source = 0.0d0 )
allocate( rtmp(mol%n*(mol%n+1)/2), source = 0.0d0 )
allocate( pbo(mol%n*(mol%n+1)/2), source = 0.0d0 )
allocate( piadr(mol%n), source = 0 )
allocate( piadr2(mol%n), source = 0 )
allocate( topo%bpair(mol%n*(mol%n+1)/2), source = 0 )
allocate( itmp(mol%n), source = 0 )
allocate( itag(mol%n), source = 0 )
allocate( sring(20,mol%n), source = 0 )
allocate( cring(10,20,mol%n), source = 0 )
allocate( piadr3(mol%n), source = 0 )
allocate( piadr4(mol%n), source = 0 )
allocate( qtmp(mol%n), source = 0.0d0 )
allocate( dxi(mol%n), source = 0.0d0 )
allocate( dgam(mol%n), source = 0.0d0 )
allocate( topo%chieeq(mol%n), source = 0.0d0 )
allocate( topo%gameeq(mol%n), source = 0.0d0 )
allocate( topo%alpeeq(mol%n), source = 0.0d0 )
allocate( topo%qa(mol%n), source = 0.0d0 )
allocate( dqa(mol%n), source = 0.0d0)
allocate( qah(mol%n), source = 0.0d0 )
allocate( nbm(20,mol%n), source = 0 )
allocate( mchar(mol%n), source = 0.0d0 )
allocate( imetal(mol%n), source = 0 )
allocate( topo%zetac6(mol%n*(mol%n+1)/2), source = 0.0d0 )
allocate( topo%xyze0(3,mol%n), source = 0.0d0 )
allocate( nbf(20,mol%n), source = 0 )
niel=0
do i=1,mol%n
niel(mol%at(i))=niel(mol%at(i))+1
enddo
write(env%unit,'(10x,"Pauling EN used:")')
do i=1,86
if(niel(i).gt.0) write(env%unit,'(10x,"Z :",i2," EN :",f6.2)') i,param%en(i)
enddo
dum = sqrt(sum(efield**2))
write(env%unit,'(10x,"electric field strengths (au):",f6.3)') dum
! alp = alp *(1.+0.0*dum)
write(env%unit,*)
write(env%unit,'(10x," ------------------------------------------------- ")')
write(env%unit,'(10x,"| Force Field Initialization |")')
write(env%unit,'(10x," ------------------------------------------------- ")')
write(env%unit,*)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! distances and bonds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
topo%xyze0 = mol%xyz ! initial geom
write(env%unit,'(10x,"distances ...")')
pbo = 0
rab = 0
sqrab = 0
do i=1,mol%n
ati=mol%at(i)
kk=i*(i-1)/2
do j=1,i-1
atj=mol%at(j)
k=kk+j
sqrab(k)=(mol%xyz(1,i)-mol%xyz(1,j))**2+(mol%xyz(2,i)-mol%xyz(2,j))**2+(mol%xyz(3,i)-mol%xyz(3,j))**2
rab(k) =sqrt(sqrab(k))
if(rab(k).lt.1.d-3) then
write(env%unit,*) i,j,ati,atj,rab(k)
call env%error("Particular close distance present", source)
exit
endif
enddo
enddo
call env%check(exitRun)
if (exitRun) then
return
end if
! Calculate CN and derivative
allocate(dcn(3,mol%n,mol%n), source = 0.0d0 )
call gfnff_dlogcoord(mol%n,mol%at,mol%xyz,rab,cn,dcn,cnthr,param) ! dcn needed
do i=1,mol%n
dum2=0
do j=1,mol%n
dum2=dum2+sqrt(dcn(1,j,i)**2+dcn(2,j,i)**2+dcn(3,j,i)**2)
enddo
mchar(i) = exp(-0.005d0*param%en(mol%at(i))**8)*dum2/(cn(i)+1.0d0) ! estimated metallic character as ratio of av. dCN and CN
! and an EN cut-off function, used in neigbor routinen and for BS estimate
enddo
deallocate(dcn)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! neighbor list, hyb and ring info
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
topo%qa = 0
qloop_count = 0
!111 continue
! do the loop only if factor is significant
do while (qloop_count.lt.2.and.gen%rqshrink.gt.1.d-3)
write(env%unit,'(10x,"----------------------------------------")')
write(env%unit,'(10x,"generating topology and atomic info file ...")')
call gfnff_neigh(env,makeneighbor,mol%n,mol%at,mol%xyz,rab,gen%rqshrink, &
& gen%rthr,gen%rthr2,gen%linthr,mchar,topo%hyb,itag,nbm,nbf,param,topo)
do i=1,mol%n
imetal(i)=param%metal(mol%at(i))
if(topo%nb(20,i).le.4.and.param%group(mol%at(i)).gt.3) imetal(i)=0 ! Sn,Pb,Bi, with small CN are better described as non-metals
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! bonds (non bonded directly in EG)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
topo%bpair=0
do i=1,mol%n
do j=1,topo%nb(20,i)
k=topo%nb(j,i)
topo%bpair(lin(k,i))=1
enddo
enddo
topo%nbond = sum(topo%bpair)
topo%nbond_blist = topo%nbond
allocate( topo%blist(2,topo%nbond), source = 0 )
allocate( btyp(topo%nbond), source = 0 )
allocate( pibo(topo%nbond), source = 0.0d0 )
pibo = -99.
topo%nbond = 0
do i=1,mol%n
kk=i*(i-1)/2
do j=1,i-1
k=kk+j
if ( topo%bpair(k) .eq. 1 ) then ! bonds
topo%nbond = topo%nbond +1
topo%blist(1,topo%nbond)=i
topo%blist(2,topo%nbond)=j
endif
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Hueckel setup for all first-row sp2 and sp atoms
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! setup list of all possible pi atoms
k=0
piadr =0
piadr2=0
do i=1,mol%n ! setup loop
piat =(topo%hyb(i).eq.1.or.topo%hyb(i).eq.2).and.pilist(mol%at(i)) ! sp or sp2 and CNOFS
kk=0
do j=1,topo%nb(20,i)
jj=topo%nb(j,i)
if(mol%at(i).eq.8.and.mol%at(jj).eq.16.and.topo%hyb(jj).eq.5) then
piat=.false.
cycle ! SO3 is not a pi
endif
if(topo%hyb(jj).eq.1.or.topo%hyb(jj).eq.2) kk=kk+1 ! attached to sp2 or sp
enddo
picon=kk.gt.0.and.nofs(mol%at(i)) ! an N,O,F (sp3) on sp2
if(mol%at(i).eq. 7.and.topo%nb(20,i).gt.3) cycle ! NR3-X is not a pi
if(mol%at(i).eq.16.and.topo%hyb(i).eq.5 ) cycle ! SO3 is not a pi
if(picon.or.piat) then
k=k+1
piadr (k)=i
piadr2(i)=k
endif
enddo
npiall=k
! make pi neighbor list
allocate( nbpi(20,npiall),pimvec(npiall), source = 0 )
nbpi=0
do i=1,mol%n
if(piadr2(i).eq.0) cycle
ii=piadr2(i)
nbpi(20,ii)=0
do j=1,topo%nb(20,i)
k=topo%nb(j,i)
if(piadr2(k).gt.0)then
nbpi(20,ii)=nbpi(20,ii)+1
nbpi(nbpi(20,ii),ii)=piadr2(k)
endif
enddo
enddo
! assign pi atoms to fragments
call mrecgff(npiall,nbpi,picount,pimvec)
deallocate(nbpi)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! setup xi correction for EEQ
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
dxi = 0 ! default none
do i=1,mol%n
ati=mol%at(i)
nn =topo%nb(20,i)
if(nn.eq.0) cycle
ip =piadr2(i)
ji =topo%nb(1,i) ! first neighbor
nh =0
nm =0
do j=1,nn
if(mol%at(topo%nb(j,i)).eq.1) nh=nh+1
if(imetal( topo%nb(j,i) ).ne.0) nm=nm+1
enddo
! hydrogen
! if(ati.eq.1.and.nn.gt.1) dxi(i)=dxi(i)-nn*0.01
! boron
if(ati.eq.5) dxi(i)=dxi(i)+nh*0.015
! carbon
if(ati.eq.6.and.nn.eq.2.and.itag(i).eq.1) dxi(i)=-0.15 ! make carbene more negative
! if(ati.eq.6.and.nn.eq.2)then
! ki=topo%nb(2,i)
! if(mol%at(ki).eq.8.and.mol%at(ji).eq.8.and.topo%nb(20,ji).eq.1.and.topo%nb(20,ki).eq.1)then ! free CO2
! dxi(ki)=0.19 ! lower EN for O
! dxi(ji)=0.19 ! " " " "
! endif
! endif
if(ati.eq.6.and.nn.eq.1.and.mol%at(ji).eq.8.and.topo%nb(20,ji).eq.1) dxi(ji)=0.15! free CO
! nitrogen
! if(ati.eq.7.and.nn.eq.1.and.mol%at(ji).eq.6) dxi(i)=0.00 !CN
! oxygen / group 6
if(ati.eq.8.and.nn.eq.1.and.ip.ne.0.and.mol%at(ji).eq.7.and.piadr2(ji).ne.0) dxi(i)= 0.05 ! nitro oxygen, otherwise NO2 HBs are too strong
if(ati.eq.8.and.nn.eq.2.and.nh.eq.2) dxi(i)=-0.02 ! H2O
if(param%group(ati).eq.6.and.nn.gt.2) dxi(i)=dxi(i)+nn*0.005! good effect
if(ati.eq.8.or.ati.eq.16) dxi(i)=dxi(i)-nh*0.005
! fluorine / group 7
if(param%group(ati).eq.7.and.ati.gt.9.and.nn.gt.1) then ! polyvalent Cl,Br ...
if(nm.eq.0)then
dxi(i)=dxi(i)-nn*0.021! good effect
else
dxi(i)=dxi(i)+nn*0.05 ! good effect for TMs
endif
endif
enddo
! prepare EEQ xi ATOMIC parameter
! at this point for the non-geom. dep. charges qa with CN = nb
do i=1,mol%n
ati=mol%at(i)
dum =min(dble(topo%nb(20,i)),gen%cnmax) ! limits it
! base val spec. corr. CN dep.
topo%chieeq(i)=-param%chi(ati) + dxi(i) + param%cnf(ati)*sqrt(dum)
topo%gameeq(i)= param%gam(ati)
if(imetal(i).eq.2)then ! the "true" charges for the TM metals are small (for various reasons)
topo%chieeq(i)=topo%chieeq(i)-gen%mchishift ! so take for the non-geom. dep. ones less electronegative metals yield more q+
endif ! which reflect better the true polarity used for guessing various
! potential terms. The positive effect of this is big.
topo%alpeeq(i)= param%alp(ati)**2
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! topology based charges
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(env%unit,'(10x,"pair mat ...")')
call nbondmat(mol%n,topo%nb,topo%bpair) ! get number of cov. bonds between atoms up to 4 bonds
write(env%unit,'(10x,"computing topology distances matrix with Floyd-Warshall algo ...")')
allocate( rabd(mol%n,mol%n), source = 0.0e0_sp)
rabd = rabd_cutoff
! determine topology distances by Floyd-Warshall algo
! they are used in the EEQ to determine qa (approximate topology charges)
do i = 1, mol%n
rabd(i, i) = 0.0
do j = 1, topo%nb(20, i)
k=topo%nb(j,i)
rabd(k, i) = param%rad(mol%at(i)) + param%rad(mol%at(k))
rabd(i, k) = rabd(k, i)
end do
end do
do k = 1, mol%n
do i = 1, mol%n
if (rabd(i, k) > gen%tdist_thr) cycle
do j = 1, mol%n
if (rabd(k, j) > gen%tdist_thr) cycle
if (rabd(i, j) > (rabd(i, k) + rabd(k, j))) then
rabd(i, j) = rabd(i, k) + rabd(k, j)
end if
end do
end do
end do
do i=1,mol%n
do j=1,i-1
ij=lin(j,i)
if(rabd(j,i).gt.gen%tdist_thr) rabd(j,i)=rabd_cutoff ! values not properly considered
rtmp(ij) = gen%rfgoed1* rabd(j,i) / 0.52917726d0
enddo
enddo
deallocate(rabd)
frag_charges_known=.false.
write(env%unit,'(10x,"making topology EEQ charges ...")')
if(topo%nfrag.le.1) then ! nothing is known
! first check for fragments
call mrecgff(mol%n,nbf,topo%nfrag,topo%fraglist)
write(env%unit,'(10x,"#fragments for EEQ constrain: ",i0)') topo%nfrag
! read QM info if it exists
call open_file(ich, 'charges', 'r')
if (ich /= -1) then
qtmp = 0
err = 0
i = 0
do while(err == 0)
read(ich,*,iostat=err) dum
if (err /= 0) exit
if (i < mol%n) then
i = i+1
qtmp(topo%fraglist(i))=qtmp(topo%fraglist(i))+dum
else
call env%warning("More charges than atoms present, assuming missmatch", source)
err = 1
end if
enddo
if (is_iostat_end(err) .and. i == mol%n) err = 0
call close_file(ich)
if (err == 0) then
if (i < mol%n .or. abs(sum(qtmp) - mol%chrg) > 1.0e-3_wp) then
call env%warning("Rejecting external charges input due to missmatch", source)
else
topo%qfrag=dnint(qtmp)
write(env%unit,'(10x,"fragment charges from <charges> :",10F7.3)') topo%qfrag(1:topo%nfrag)
end if
else
call env%warning("Could not initialize fragment charges from file", source)
endif
call env%check(exitRun)
if (exitRun) then
return
end if
endif
if(mol%n.lt.100.and.topo%nfrag.gt.2.and.mol%chrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
itmp=0
do i=1,mol%n
itmp(topo%fraglist(i))=itmp(topo%fraglist(i))+1
enddo
do i=1,topo%nfrag
write(env%unit,*)i,itmp(i)
enddo
call env%error('fragment charge input required', source)
return
endif
if(mol%n.ge.100.and.topo%nfrag.gt.2.and.mol%chrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
topo%qfrag(1)=mol%chrg
topo%qfrag(2:topo%nfrag)=0
endif
if(topo%nfrag.eq.2.and.mol%chrg.ne.0.and.sum(topo%qfrag(2:topo%nfrag)).gt.999) then
write(env%unit,*) 'trying auto detection of charge on 2 fragments:'
topo%qfrag(1)=0
topo%qfrag(2)=mol%chrg
call goedeckera(env,mol%n,mol%at,topo%nb,rtmp,topo%qa,dum1,topo)
call env%check(exitRun)
if (exitRun) then
call env%error("Failed to generate charges", source)
return
end if
topo%qfrag(2)=0
topo%qfrag(1)=mol%chrg
call goedeckera(env,mol%n,mol%at,topo%nb,rtmp,topo%qa,dum2,topo)
call env%check(exitRun)
if (exitRun) then
call env%error("Failed to generate charges", source)
return
end if
if(dum1.lt.dum2) then
topo%qfrag(1)=0
topo%qfrag(2)=mol%chrg
endif
write(env%unit,*) 'dEes :',dum1-dum2
write(env%unit,*) 'charge 1/2:',topo%qfrag(1:2)
endif
else if (allocated(mol%pdb)) then ! frag_charges_known
write(env%unit,'(10x,"#fragments for EEQ constrain from pdb file: ",i0)') topo%nfrag
frag_charges_known=.true.
endif
! make estimated, topology only EEQ charges from rabd values, including "right" fragment charge
call goedeckera(env,mol%n,mol%at,topo%nb,rtmp,topo%qa,ees,topo)
call env%check(exitRun)
if (exitRun) then
call env%error("Failed to generate charges", source)
return
end if
! estimate how much of the frag charge is on the pi sub systems
if(picount.gt.0.and.qloop_count.gt.0) then
allocate( ipis(picount), source = 0 )
if(frag_charges_known) then ! PDB case
ipis = 0
do pis=1,picount ! loop over pi systems
do k=1,npiall
if(pimvec(k).eq.pis) then
ipis(pis)=ipis(pis)+topo%qpdb(piadr(k))
topo%qpdb(piadr(k))=0
endif
enddo
enddo
else ! general case
qtmp = topo%qa ! save the "right" ones
qah = topo%qa
call qheavy(mol%n,mol%at,topo%nb,qah) ! heavy atoms only ie H condensed to neighbor
do pis=1,picount ! loop over pi systems
do k=1,npiall
if(pimvec(k).eq.pis) then
kk=piadr(k)
ifrag=topo%fraglist(kk) !the pi atom of this pi fragment is in EEQ fragment ifrag
exit
endif
enddo
dum2=topo%qfrag(ifrag) ! save
topo%qfrag(ifrag) = 0 ! make only this EEQ fragment neutral
call goedeckera(env,mol%n,mol%at,topo%nb,rtmp,topo%qa,ees,topo) ! for neutral
call env%check(exitRun)
if (exitRun) then
call env%error("Failed to generate charges", source)
return
end if
topo%qfrag(ifrag) = dum2 ! back
call qheavy(mol%n,mol%at,topo%nb,topo%qa)
dqa =qah-topo%qa ! difference charges upon ionization
dum1=0
dum=0
do k=1,npiall
if(pimvec(k).eq.pis) dum=dum+dqa(piadr(k)) ! only pi atoms
enddo
dum = dum * 1.1 !charges tend to be slightly too small 1.1-1.2
ipis(pis)=idnint(dum)
dum1=dum1+dum
enddo
topo%qa = qtmp ! put "right" charges used in FF construction and for HB/XB in place
endif
endif
if(qloop_count.eq.0) itmp(1:mol%n)=topo%nb(20,1:mol%n)
qloop_count=qloop_count+1
if(qloop_count.lt.2.and.gen%rqshrink.gt.1.d-3) then ! do the loop only if factor is significant
deallocate(topo%blist,btyp,pibo,pimvec)
! goto 111
endif
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! change EEQ J with estimated q
! which is a kind of third-order term
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,mol%n
ff=0 ! do nothing
if(mol%at(i).eq. 1) ff=-0.08 ! H
if(mol%at(i).eq. 5) ff=-0.05 ! B
if(mol%at(i).eq. 6) then
ff=-0.27 ! C
if(topo%hyb(i).lt.3) ff=-0.45 ! unsat
if(topo%hyb(i).lt.2) ff=-0.34 ! unsat
endif
if(mol%at(i).eq. 7) then
ff=-0.13 ! N
if(piadr(i).ne.0) ff=-0.14
if(amide(mol%n,mol%at,topo%hyb,topo%nb,piadr,i))ff=-0.16
endif
if(mol%at(i).eq. 8) then
ff=-0.15 ! O
if(topo%hyb(i).lt.3) ff=-0.08 ! unsat
endif
if(mol%at(i).eq. 9) ff= 0.10 ! F
if(mol%at(i).gt.10) ff=-0.02 ! heavy
if(mol%at(i).eq.17) ff=-0.02 ! Cl
if(mol%at(i).eq.35) ff=-0.11 ! Br
if(mol%at(i).eq.53) ff=-0.07 ! I
if(imetal(i).eq.1) ff=-0.08 ! M maing
if(imetal(i).eq.2) ff=-0.9 ! M TM ??? too large
if(param%group(mol%at(i)).eq.8) ff= 0.0 ! RG
dgam(i)=topo%qa(i)*ff
enddo
! prepare true EEQ parameter, they are ATOMIC not element specific!
do i=1,mol%n
! base val spec. corr.
topo%chieeq(i)=-param%chi(mol%at(i)) + dxi(i)
topo%gameeq(i)= param%gam(mol%at(i)) +dgam(i)
if (amideH(mol%n,mol%at,topo%hyb,topo%nb,piadr2,i)) topo%chieeq(i) = topo%chieeq(i) - 0.02
ff = 0
if(mol%at(i).eq.6) ff= 0.09
if(mol%at(i).eq.7) ff=-0.21
if(param%group(mol%at(i)).eq.6)ff=-0.03
if(param%group(mol%at(i)).eq.7)ff= 0.50
if(imetal(i).eq.1) ff= 0.3
if(imetal(i).eq.2) ff=-0.1
topo%alpeeq(i) = (param%alp(mol%at(i))+ff*topo%qa(i))**2
enddo
deallocate(dgam,dxi)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! get ring info (smallest ring size)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(env%unit,'(10x,"rings ...")')
!$omp parallel default(none) private(i,cr,sr) shared(mol,nbm,cring,sring)
!$omp do
do i=1,mol%n
call getring36(mol%n,mol%at,nbm,i,cr,sr)
cring(1:10,1:20,i)=cr(1:10,1:20)
sring( 1:20,i)=sr(1:20)
enddo
!$omp end do
!$omp end parallel
deallocate(nbm)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! bonded atom triples not included in
! bend and tors
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
idum=1000*mol%n
allocate( topo%b3list(3,idum), source = 0 )
topo%nbatm=0
do i=1,mol%n
do j=1,i-1
ij=lin(j,i)
if(topo%bpair(ij).eq.3) then ! 1,4 exclusion of back-pair makes it worse, 1,3 makes little effect
do m=1,topo%nb(20,j)
k=topo%nb(m,j)
topo%nbatm=topo%nbatm+1
topo%b3list(1,topo%nbatm)=i
topo%b3list(2,topo%nbatm)=j
topo%b3list(3,topo%nbatm)=k
enddo
do m=1,topo%nb(20,i)
k=topo%nb(m,i)
topo%nbatm=topo%nbatm+1
topo%b3list(1,topo%nbatm)=i
topo%b3list(2,topo%nbatm)=j
topo%b3list(3,topo%nbatm)=k
enddo
endif
enddo
enddo
if(topo%nbatm.gt.idum) then
write(env%unit,*) idum,topo%nbatm
call env%error('overflow in ini', source)
return
endif
write(env%unit,'(10x,"# BATM",3x,i0)') topo%nbatm
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! non bonded pair exponents
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,mol%n
ati=mol%at(i)
fn=1.0d0 + gen%nrepscal/(1.0d0+dble(topo%nb(20,i))**2)
dum1=param%repan(ati)*(1.d0 + topo%qa(i)*gen%qrepscal)*fn ! a small but physically correct decrease of repulsion with q
f1=zeta(ati,topo%qa(i))
do j=1,i-1
atj=mol%at(j)
fn=1.0d0 + gen%nrepscal/(1.0d0+dble(topo%nb(20,j))**2)
dum2=param%repan(atj)*(1.d0 + topo%qa(j)*gen%qrepscal)*fn
f2=zeta(atj,topo%qa(j))
ij=lin(j,i)
ff = 1.0d0
if(ati.eq.1.and.atj.eq.1) then
ff = 1.0d0*gen%hhfac ! special H ... H case (for other pairs there is no good effect of this)
if(topo%bpair(ij).eq.3) ff=ff*gen%hh14rep ! 1,4 case important for right torsion pot.
if(topo%bpair(ij).eq.2) ff=ff*gen%hh13rep ! 1,3 case
endif
if((ati.eq.1.and.param%metal(atj).gt.0).or.(atj.eq.1.and.param%metal(ati).gt.0)) ff=0.85 ! M...H
if((ati.eq.1.and.atj.eq.6).or.(atj.eq.1.and.ati.eq.6)) ff=0.91 ! C...H, good effect
if((ati.eq.1.and.atj.eq.8).or.(atj.eq.1.and.ati.eq.8)) ff=1.04 ! O...H, good effect
topo%alphanb(ij)=sqrt(dum1*dum2)*ff
topo%zetac6(ij)=f1*f2 ! D4 zeta scaling using qref=0
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! make list of HB donor bascity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!atom specific (not element) basicity parameters
allocate( topo%hbbas(mol%n), source =1.0d0 )
do i = 1,mol%n
nn=topo%nb(20,i)
ati=mol%at(i)
topo%hbbas(i)=param%xhbas(mol%at(i))
! Carbene:
if(ati.eq.6.and.nn.eq.2.and.itag(i).eq.1) topo%hbbas(i) = 1.46
! Carbonyl R-C=O
if(ati.eq.8.and.nn.eq.1.and.mol%at(topo%nb(nn,i)).eq.6) topo%hbbas(i) = 0.68
! Nitro R-N=O
if(ati.eq.8.and.nn.eq.1.and.mol%at(topo%nb(nn,i)).eq.7) topo%hbbas(i) = 0.47
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! make list of HB donor acidity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!atom specific (not element) basicity parameters
allocate( topo%hbaci(mol%n), source =1.0d0 )
do i = 1,mol%n
topo%hbaci(i)=param%xhaci(mol%at(i))
end do
do i = 1,mol%n
nn=topo%nb(1,i)
topo%hbaci(i)=param%xhaci(mol%at(i))
! AmideH:
if (amideH(mol%n,mol%at,topo%hyb,topo%nb,piadr2,i)) topo%hbaci(nn) = topo%hbaci(nn) * 0.80
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! make list of ABs for HAB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
allocate( topo%hbatHl(mol%n),topo%hbatABl(2,mol%n*(mol%n+1)/2), source = 0 )
topo%nathbH=0
do i=1,mol%n
if(mol%at(i).ne.1) cycle
if(topo%hyb(i).eq.1) cycle ! exclude bridging hydrogens from HB correction
ff=gen%hqabthr
j=topo%nb(1,i)
if(j.le.0) cycle
if(mol%at(j).gt.10) ff=ff-0.20 ! H on heavy atoms may be negatively charged
if(mol%at(j).eq.6.and.topo%hyb(j).eq.3) ff=ff+0.05 ! H on sp3 C must be really positive 0.05
if(topo%qa(i).gt.ff)then ! make list of HB H atoms but only if they have a positive charge
topo%nathbH=topo%nathbH+1
topo%hbatHl(topo%nathbH)=i
endif
enddo
write(env%unit,'(10x,"# H in HB",3x,i0)') topo%nathbH
topo%nathbAB=0
do i=1,mol%n
if(mol%at(i).eq. 6.and.piadr2(i).eq.0) cycle ! C sp or sp2 pi
ff=gen%qabthr
if(mol%at(i).gt.10) ff=ff+0.2 ! heavy atoms may be positively charged
if(topo%qa(i).gt.ff) cycle
do j=1,i-1
ff=gen%qabthr
if(mol%at(j).gt.10) ff=ff+0.2 ! heavy atoms may be positively charged
if(topo%qa(j).gt.ff) cycle
call hbonds(i,j,hbpi,hbpj,param,topo)
if(hbpi(1)*hbpj(2).lt.1.d-6.and.hbpi(2)*hbpj(1).lt.1.d-6)cycle
if(mol%at(j).eq. 6.and.piadr2(j).eq.0) cycle ! C sp or sp2 pi
topo%nathbAB = topo%nathbAB + 1
topo%hbatABl(1,topo%nathbAB)=i
topo%hbatABl(2,topo%nathbAB)=j
enddo
enddo
! make ABX list
m=0
do i=1,mol%n
do ia=1,topo%nb(20,i)
ix=topo%nb(ia,i)
if(xatom(mol%at(ix))) then
if(mol%at(ix).eq.16.and.topo%nb(20,ix).gt.2) cycle ! no sulphoxide etc S
do j=1,mol%n
if(i.eq.j.or.j.eq.ix) cycle
if(topo%bpair(lin(j,ix)).le.3) cycle ! must be A...B and not X-B i.e. A-X...B
if(param%xhbas(mol%at(j)).lt.1.d-6) cycle ! B must be O,N,...
if(param%group(mol%at(j)).eq.4 ) then
if(piadr2(j).eq.0.or.topo%qa(j).gt.0.05) cycle ! must be a (pi)base
endif
m=m+1
enddo
endif
enddo
enddo
topo%natxbAB=m
allocate(topo%xbatABl(3,topo%natxbAB), source = 0 )
m=0
do i=1,mol%n
do ia=1,topo%nb(20,i)
ix=topo%nb(ia,i)
if(xatom(mol%at(ix))) then
if(mol%at(ix).eq.16.and.topo%nb(20,ix).gt.2) cycle ! no sulphoxide etc S
do j=1,mol%n
if(i.eq.j.or.j.eq.ix) cycle
if(topo%bpair(lin(j,ix)).le.3) cycle ! must be A...B and not X-B i.e. A-X...B
if(param%xhbas(mol%at(j)).lt.1.d-6) cycle ! B must be O,N,...
if(param%group(mol%at(j)).eq.4 ) then
if(piadr2(j).eq.0.or.topo%qa(j).gt.0.05) cycle ! must be a (pi)base
endif
m=m+1
topo%xbatABl(1,m)=i
topo%xbatABl(2,m)=j
topo%xbatABl(3,m)=ix
enddo
endif
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! do Hueckel
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(picount.gt.0) then
write(env%unit,'(10x,"doing iterative Hueckel for ",i0," subsystem(s) ...")') picount
allocate( pispop(picount),pisip(picount),pisea(picount), source = 0.0d0 )
allocate( piel(mol%n), source = 0 )
itmp = 0 ! save pi atom info
hcalc= 0
pisip= 0
pisea= 0
if(pr) then
write(env%unit,'(10x,"iterative Hueckel run to get P ...")')
endif
do pis=1,picount ! loop over pi systems
npi =0
nelpi =0
piadr3=0
piadr4=0
piel =0
do k=1,npiall
if(pimvec(k).eq.pis) then
npi=npi+1
ati =mol%at (piadr(k))
hybi=topo%hyb(piadr(k))
ii =nelpi
if(ati.eq.5.and.hybi.eq.1) nelpi=nelpi+1 ! B in borine
if(ati.eq.6.and.itag(piadr(k)).ne.1) nelpi=nelpi+1 ! skip if its a carbene (tag itag=1)
if(ati.eq.7.and.hybi.eq.2.and.itag(piadr(k)).eq.1) &
& nelpi=nelpi+1 ! the itag=1 avoids an odd el number for the nitro group (its 4)
if(ati.eq.7.and.hybi.le.2) nelpi=nelpi+1
if(ati.eq.7.and.hybi.eq.3) nelpi=nelpi+2
if(ati.eq.8.and.hybi.eq.1) nelpi=nelpi+1
if(ati.eq.8.and.hybi.eq.2) nelpi=nelpi+1
if(ati.eq.8.and.hybi.eq.3) nelpi=nelpi+2
if(ati.eq.9.and.hybi.ne.1) nelpi=nelpi+2
if(ati.eq.9.and.hybi.eq.1) nelpi=nelpi+3 !??? otherwise fluor-furan+ is wrong
if(ati.eq.16.and.hybi.eq.1) nelpi=nelpi+1
if(ati.eq.16.and.hybi.eq.2) nelpi=nelpi+1
if(ati.eq.16.and.hybi.eq.3) nelpi=nelpi+2
if(ati.eq.17.and.hybi.eq.0) nelpi=nelpi+2
if(ati.eq.17.and.hybi.eq.1) nelpi=nelpi+3
piadr3(npi)=piadr(k) ! map to original, full atom set
piadr4(piadr(k))=npi
piel(piadr(k))=nelpi-ii
if(piel(piadr(k)).gt.2)piel(piadr(k))=2
endif
enddo
nelpi=nelpi-ipis(pis)
if(npi.lt.2.or.nelpi.lt.1) cycle
allocate(Api(npi,npi),apisave(npi,npi),Pold(npi,npi),S(npi,npi),occ(npi),eps(npi)) ! S is just scratch here
eold= 0
Pold= 2.d0/3.d0
! iterative Hueckel loop, off-diag terms are reduced depending on P to avoid overdelocalization
do nn=1,nint(gen%maxhiter) ! just some iterations
Api = 0
do i=1,npi
ii=piadr3(i)
Api(i,i)=gen%hdiag(mol%at(ii))+topo%qa(ii)*gen%hueckelp3-dble(piel(ii)-1)*gen%pilpf
enddo
! loop over bonds for pair interactions
do i=1,topo%nbond
ii=topo%blist(1,i)
jj=topo%blist(2,i)
ia=piadr4(ii)
ja=piadr4(jj)
if(ia.gt.0.and.ja.gt.0) then
dum=1.d-9*rab(lin(ii,jj)) ! distort so that Huckel for e.g. COT localizes to right bonds
dum=sqrt(gen%hoffdiag(mol%at(ii))*gen%hoffdiag(mol%at(jj)))-dum ! better than arithmetic
dum2=gen%hiter
if(topo%hyb(ii).eq.1) dum2=dum2*gen%htriple ! triple bond is different
if(topo%hyb(jj).eq.1) dum2=dum2*gen%htriple ! triple bond is different
Api(ja,ia)=-dum * (1.0d0-dum2*(2.0d0/3.0d0-Pold(ja,ia))) ! Pmat scaling with benzene as reference
Api(ia,ja)=Api(ja,ia)
endif
enddo
apisave = Api
call gfnffqmsolve(.false.,Api,S,.false.,4000.0d0,npi,0,nelpi,dum,occ,eps) !diagonalize, 4000 better than 300
do i=1,npi ! save IP/EA
if(occ(i).gt.0.5) then
pisip(pis)=eps(i) ! IP
if(i+1.lt.npi)pisea(pis)=eps(i+1) ! EA
endif
enddo
if(abs(dum-eold).lt.1.d-4) exit ! end of iterations
Pold = Api
eold = dum
enddo
! end of iterative loop
if(pr)then
write(env%unit,'(''Hueckel system :'',i3,'' charge : '',i3,'' ndim/Nel :'',2i5, &
& 3x, ''eps(HOMO/LUMO)'',2f12.6)')pis,ipis(pis),npi,nelpi,pisip(pis),pisea(pis)
end if
if(pisip(pis).gt.0.40) then
write(env%unit,'(a,i0,a)')'WARNING: probably wrong pi occupation for system ',pis,'. Second attempt with Nel=Nel-1!'
do i=1,mol%n
if(piadr4(i).ne.0) write(env%unit,*) 'at,nb,topo%hyb,Npiel:', i,mol%sym(i),topo%nb(20,i),topo%hyb(i),piel(i)
enddo
nelpi=nelpi-1
Api = Apisave
call gfnffqmsolve(.false.,Api,S,.false.,4000.0d0,npi,0,nelpi,dum,occ,eps) !diagonalize
call PREIG(6,occ,1.0d0,eps,1,npi)
do i=1,npi ! save IP/EA
if(occ(i).gt.0.5) then
pisip(pis)=eps(i) ! IP
if(i+1.lt.npi)pisea(pis)=eps(i+1) ! EA
endif
enddo
if(pr)then
write(env%unit,'(''Hueckel system :'',i3,'' charge : '',i3,'' ndim/Nel :'',2i5, &
& 3x, ''eps(HOMO/LUMO)'',2f12.6)')pis,ipis(pis),npi,nelpi,pisip(pis),pisea(pis)
end if
endif
! save BO
do i=1,topo%nbond
ii=topo%blist(1,i)
jj=topo%blist(2,i)
ia=piadr4(ii)
ja=piadr4(jj)
if(ia.gt.0.and.ja.gt.0) then
pibo(i)=Api(ja,ia)
pbo(lin(ii,jj))=Api(ja,ia)
itmp(ii)=1
itmp(jj)=1
endif
enddo
deallocate(Api,apisave,Pold,S,occ,eps)
enddo
! end of pi system loop
piadr = itmp ! array used for identifying pi atoms in following codes
deallocate(pispop,pisip,pisea,ipis,pimvec,piel)
endif
!----------- end Hueckel
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! modify hyb due to pi assignment
! and output
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=1,mol%n
if(topo%hyb(i).eq.2 .and. piadr(i).eq.0 .and. topo%nb(20,i).eq.3 .and. param%group(mol%at(i)).eq.4)then ! C,Si,Ge... CN=3, no pi
jj=topo%nb(1,i)
kk=topo%nb(2,i)
ll=topo%nb(3,i)
phi=omega(mol%n,mol%xyz,i,jj,kk,ll) ! the shitty second geom. dep. term GEODEP
if(abs(phi)*180./pi.gt.40.d0) topo%hyb(i) = 3 ! change to sp^3
endif
enddo
! if(pr)then
write(env%unit,*)
write(env%unit,'(2x,"atom neighbors erfCN metchar sp-hybrid imet pi qest coordinates")')
do i=1,mol%n
j = topo%hyb(i)
if(amide(mol%n,mol%at,topo%hyb,topo%nb,piadr,i)) j=-topo%hyb(i)
if(mol%at(i).eq.6.and.itag(i).eq.1) j=-topo%hyb(i)
write(env%unit,'(i5,2x,a2,3x,i4,3x,f5.2,2x,f5.2,8x,i2,3x,i2,3x,i2,2x,f6.3,3f12.6)') &
& i,mol%sym(i),topo%nb(20,i),cn(i),mchar(i),j,imetal(i),piadr(i),topo%qa(i),mol%xyz(1:3,i)
enddo
! compute fragments and charges for output (check for CT)
! call mrecgff(mol%n,topo%nb,nmol,piadr3)
! write(env%unit,*) 'Nmol',nmol
if(pr)then
write(env%unit,'(/,''molecular fragment # atoms topo charge'')')
do i=1,topo%nfrag
dum=0
m=0
do k=1,mol%n
if(topo%fraglist(k).eq.i) then
m=m+1
dum=dum+topo%qa(k)