1 subroutine w2x(nvarr,x,*)
4 include "DIMENSIONS.ZSCOPT"
5 include "COMMON.CONTROL"
6 include "COMMON.INTERACT"
7 include "COMMON.TORSION"
9 include "COMMON.WEIGHTS"
10 include "COMMON.ENERGIES"
11 include "COMMON.IOUNITS"
12 include "COMMON.NAMES"
14 include "COMMON.VMCPAR"
15 include "COMMON.CLASSES"
16 include "COMMON.WEIGHTDER"
17 include "COMMON.FFIELD"
18 include "COMMON.SCCOR"
19 integer itor2typ(-2:2) /-20,-20,10,9,20/
20 integer ind_shield /25/
22 double precision x(max_paropt),xchuj,epsi
23 integer nvarr,i,ii,j,k,l,ll,iprot,iblock
28 if (imask(i).eq.1 .and. i.ne.ind_shield) then
37 c Change 9/21/07 AL Variables are epsilons and not the increments on epsilons
38 c to enable control on their deviations from original values.
39 c Single-body epsilons are now defined as epsilons of pairs of same type side
42 c write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc,
43 c & " npair_sc",npair_sc
48 eps_orig(i,j)=eps(i,j)
56 c write (iout,*) "ii",ii," x",x(ii)
59 c if (epss_up(i).eq.0.0d0) then
60 c x_up(ii)=dabs(epsi)*epss_low(i)
61 c x_low(ii)=epsi-x_up(ii)
62 c x_up(ii)=epsi+x_up(ii)
64 c x_low(ii)=epss_low(i)
74 if (ityp_ssc(j).eq.iti) then
75 write (iout,*) "Error - pair ",restyp(iti),"-",
77 & " specified in variable epsilon list but type defined",
78 & " in single-type list."
87 c if (epsp_up(i).eq.0.0d0) then
88 c x_up(ii)=dabs(epsi)*epsp_low(i)
89 c x_low(ii)=epsi-x_up(ii)
90 c x_up(ii)=epsi+x_up(ii)
92 c x_low(ii)=epsp_low(i)
95 c write (iout,*) "ii",ii," x",x(ii)
105 write (iout,*) "Indices of optimizable torsionals"
106 do i=-ntortyp+1,ntortyp-1
107 do j=-ntortyp+1,ntortyp-1
108 c write (iout,*) "mask",i,j,mask_tor(j,i)
110 if (tor_mode.eq.0) then
113 if (mask_tor(0,j,i,iblock).eq.1) then
116 x(ii)=weitor(0,j,i,iblock)
117 x_low(ii)=weitor_low(0,j,i,iblock)
118 x_up(ii)=weitor_up(0,j,i,iblock)
119 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
120 & restyp(itor2typ(i)),iblock,ii
121 else if (mask_tor(0,j,i,iblock).eq.2) then
122 do k=1,nterm(j,i,iblock)
125 x(ii)=v1(k,j,i,iblock)
126 x_low(ii)=weitor_low(0,j,i,iblock)
127 x_up(ii)=weitor_up(0,j,i,iblock)
128 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
129 & restyp(itor2typ(i)),iblock,k,ii
131 else if (mask_tor(0,j,i,iblock).gt.2) then
132 do k=1,nterm(j,i,iblock)
135 x(ii)=v1(k,j,i,iblock)
136 x_low(ii)=weitor_low(0,j,i,iblock)
137 x_up(ii)=weitor_up(0,j,i,iblock)
138 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
139 & restyp(itor2typ(i)),iblock,k,ii
141 do k=1,nterm(j,i,iblock)
144 x(ii)=v2(k,j,i,iblock)
145 x_low(ii)=weitor_low(0,j,i,iblock)
146 x_up(ii)=weitor_up(0,j,i,iblock)
147 write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
148 & restyp(itor2typ(i)),ii,k,ii
153 else if (tor_mode.eq.1) then
155 if (mask_tor(0,j,i,1).eq.1) then
158 x(ii)=weitor(0,j,i,1)
159 x_low(ii)=weitor_low(0,j,i,1)
160 x_up(ii)=weitor_up(0,j,i,1)
161 write (iout,'(2a4,i4)') restyp(itor2typ(j)),
162 & restyp(itor2typ(i)),ii
163 else if (mask_tor(0,j,i,1).eq.2) then
164 do k=1,nterm_kcc(j,i)
165 do l=1,nterm_kcc_Tb(j,i)
166 do ll=1,nterm_kcc_Tb(j,i)
169 x(ii)=v1_kcc(ll,l,k,j,i)
170 x_low(ii)=weitor_low(0,j,i,1)
171 x_up(ii)=weitor_up(0,j,i,1)
172 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
173 & restyp(itor2typ(i)),k,ii
177 else if (mask_tor(0,j,i,1).gt.2) then
178 do k=1,nterm_kcc(j,i)
179 do l=1,nterm_kcc_Tb(j,i)
180 do ll=1,nterm_kcc_Tb(j,i)
183 x(ii)=v1_kcc(ll,l,k,j,i)
184 x_low(ii)=weitor_low(0,j,i,1)
185 x_up(ii)=weitor_up(0,j,i,1)
186 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
187 & restyp(itor2typ(i)),k,ii
191 do k=1,nterm_kcc(j,i)
192 do l=1,nterm_kcc_Tb(j,i)
193 do ll=1,nterm_kcc_Tb(j,i)
196 x(ii)=v2_kcc(ll,l,k,j,i)
197 x_low(ii)=weitor_low(0,j,i,1)
198 x_up(ii)=weitor_up(0,j,i,1)
199 write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
200 & restyp(itor2typ(i)),k,ii
207 c Compute torsional coefficients from local energy surface expansion coeffs
208 if (mask_tor(0,j,i,1).eq.1) then
210 write (iout,*) "i",i," j",j," weitor",weitor(0,j,i,1)
211 x(ii)=weitor(0,j,i,1)
212 x_low(ii)=weitor_low(0,j,i,1)
213 x_up(ii)=weitor_up(0,j,i,1)
216 write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
225 nbacktor_var=ntor_var
226 write (iout,*) "nbacktor_var",nbacktor_var
230 write (iout,*) "Indices of optimizable SCCOR torsionals"
231 do i=-nsccortyp,nsccortyp
232 do j=-nsccortyp,nsccortyp
234 c write (iout,*) "mask",i,j,mask_tor(j,i)
235 if (mask_tor(k,j,i,1).eq.1) then
238 x(ii)=weitor(k,j,i,1)
239 x_low(ii)=weitor_low(k,j,i,1)
240 x_up(ii)=weitor_up(k,j,i,1)
241 write (iout,'(i4,2a4,i4)') k,restyp(j),restyp(i),ii
242 else if (mask_tor(k,j,i,1).eq.2) then
243 do l=1,nterm_sccor(j,i)
246 x(ii)=v1sccor(l,k,j,i)
247 x_low(ii)=weitor_low(k,j,i,1)
248 x_up(ii)=weitor_up(k,j,i,1)
249 write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
251 else if (mask_tor(k,j,i,1).gt.2) then
252 do l=1,nterm_sccor(j,i)
255 x(ii)=v1sccor(l,k,j,i)
256 x_low(ii)=weitor_low(k,j,i,1)
257 x_up(ii)=weitor_up(k,j,i,1)
258 write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
260 do l=1,nterm_sccor(j,i)
263 x(ii)=v2sccor(l,k,j,i)
264 x_low(ii)=weitor_low(k,j,i,1)
265 x_up(ii)=weitor_up(k,j,i,1)
266 write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
275 c 7/8/17 AL Optimization of the bending potentials
280 if (mask_ang(i).eq.0) cycle
281 do j=1,nbend_kcc_TB(i)
284 x(ii)=v1bend_chyb(j,i)
285 x_low(ii)=v1bend_low(j,i)
286 x_up(ii)=v1bend_up(j,i)
290 write (iout,*) "nang_var",nang_var
292 c Optimized cumulant coefficients
294 if (mod_fourier(nloctyp).gt.0) then
299 if (mod_fourier(i).gt.0) then
303 if (mask_bnew1(k,j,i).gt.0) then
306 x_low(ii)=bnew1_low(k,j,i)
307 x_up(ii)=bnew1_up(k,j,i)
314 if (mask_bnew2(k,j,i).gt.0) then
317 x_low(ii)=bnew2_low(k,j,i)
318 x_up(ii)=bnew2_up(k,j,i)
325 if (mask_ccnew(k,j,i).gt.0) then
328 x_low(ii)=ccnew_low(k,j,i)
329 x_up(ii)=ccnew_up(k,j,i)
336 if (mask_ddnew(k,j,i).gt.0) then
339 x_low(ii)=ddnew_low(k,j,i)
340 x_up(ii)=ddnew_up(k,j,i)
346 if (mask_e0new(j,i).gt.0) then
349 x_low(ii)=e0new_low(j,i)
350 x_up(ii)=e0new_up(j,i)
357 if (mask_eenew(l,k,j,i).gt.0) then
360 x_low(ii)=eenew_low(l,k,j,i)
361 x_up(ii)=eenew_up(l,k,j,i)
371 IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN
375 if (mod_fouriertor(i).gt.0) then
379 if (mask_bnew1tor(k,j,i).gt.0) then
381 x(ii)=bnew1tor(k,j,i)
382 x_low(ii)=bnew1tor_low(k,j,i)
383 x_up(ii)=bnew1tor_up(k,j,i)
390 if (mask_bnew2tor(k,j,i).gt.0) then
392 x(ii)=bnew2tor(k,j,i)
393 x_low(ii)=bnew2tor_low(k,j,i)
394 x_up(ii)=bnew2tor_up(k,j,i)
401 if (mask_ccnewtor(k,j,i).gt.0) then
403 x(ii)=ccnewtor(k,j,i)
404 x_low(ii)=ccnewtor_low(k,j,i)
405 x_up(ii)=ccnewtor_up(k,j,i)
412 if (mask_ddnewtor(k,j,i).gt.0) then
414 x(ii)=ddnewtor(k,j,i)
415 x_low(ii)=ddnewtor_low(k,j,i)
416 x_up(ii)=ddnewtor_up(k,j,i)
422 if (mask_e0newtor(j,i).gt.0) then
425 x_low(ii)=e0newtor_low(j,i)
426 x_up(ii)=e0newtor_up(j,i)
433 if (mask_eenewtor(l,k,j,i).gt.0) then
435 x(ii)=eenewtor(l,k,j,i)
436 x_low(ii)=eenewtor_low(l,k,j,i)
437 x_up(ii)=eenewtor_up(l,k,j,i)
451 if (mod_fourier(i).gt.0) then
453 if (mask_fourier(j,i).gt.0) then
466 c Added SC sigmas and anisotropies
468 if (mod_side_other) then
470 if (mask_sigma(1,i).gt.0) then
473 x_low(ii)=sigma_low(1,i)
474 x_up(ii)=sigma_up(1,i)
475 c write (iout,*) "sig0",i,ii,x(ii), x_low(ii),x_up(ii)
477 if (mask_sigma(2,i).gt.0) then
480 x_low(ii)=sigma_low(2,i)
481 x_up(ii)=sigma_up(2,i)
482 c write (iout,*) "sigi",i,ii,x(ii), x_low(ii),x_up(ii)
484 if (mask_sigma(3,i).gt.0) then
487 x_low(ii)=sigma_low(3,i)
488 x_up(ii)=sigma_up(3,i)
489 c write (iout,*) "chi",i,ii,x(ii), x_low(ii),x_up(ii)
491 if (mask_sigma(4,i).gt.0) then
494 x_low(ii)=sigma_low(4,i)
495 x_up(ii)=sigma_up(4,i)
496 c write (iout,*) "alp",i,ii,x(ii), x_low(ii),x_up(ii)
506 if (mask_elec(i,j,1).gt.0) then
509 x_low(ii)=epp_low(i,j)
512 if (mask_elec(i,j,2).gt.0) then
515 x_low(ii)=rpp_low(i,j)
518 if (mask_elec(i,j,3).gt.0) then
521 x_low(ii)=elpp6_low(i,j)
522 x_up(ii)=elpp6_up(i,j)
524 if (mask_elec(i,j,4).gt.0) then
527 x_low(ii)=elpp3_low(i,j)
528 x_up(ii)=elpp3_up(i,j)
535 C Define the SC-p interaction constants
539 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
540 & mask_scp(0,2,2) .gt. 0) then
542 if (mask_scp(0,j,1).gt.0) then
545 x_low(ii)=epscp_low(0,j)
546 x_up(ii)=epscp_up(0,j)
548 if (mask_scp(0,j,2).gt.0) then
551 x_low(ii)=rscp_low(0,j)
552 x_up(ii)=rscp_up(0,j)
558 if (mask_scp(i,j,1).gt.0) then
561 x_low(ii)=epscp_low(i,j)
562 x_up(ii)=epscp_up(i,j)
564 if (mask_scp(i,j,2).gt.0) then
567 x_low(ii)=rscp_low(i,j)
568 x_up(ii)=rscp_up(i,j)
576 c Weight of the shielding term
577 if (imask(ind_shield).eq.1) then
580 x_low(ii)=ww_low(ind_shield)
581 x_up(ii)=ww_up(ind_shield)
585 c Compute boundaries, if defined relative to starting values
587 write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i)
588 if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then
592 else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then
594 if (x(i).gt.0.0d0) then
595 x_low(i)=x(i)*(1.0d0-xchuj)
596 x_up(i)=x(i)*(1.0d0+xchuj)
598 x_low(i)=x(i)*(1.0d0+xchuj)
599 x_up(i)=x(i)*(1.0d0-xchuj)
604 c Base of heat capacity
607 if (nbeta(2,iprot).gt.0) then
609 x(ii)=heatbase(iprot)
617 write (iout,*) "Number of optimizable parameters:",nvarr
619 write (iout,*) "Initial variables and boundaries:"
620 write (iout,'(3x,3a15)') " x"," x_low"," x_up"
622 write (iout,'(i3,3e15.5)') i,x(i),x_low(i),x_up(i)
627 c------------------------------------------------------------------------------
628 subroutine x2w(nvarr,x)
631 include "DIMENSIONS.ZSCOPT"
632 include "COMMON.CONTROL"
633 include "COMMON.NAMES"
634 include "COMMON.WEIGHTS"
635 include "COMMON.INTERACT"
636 include "COMMON.ENERGIES"
637 include "COMMON.FFIELD"
638 include "COMMON.TORSION"
639 include "COMMON.LOCAL"
640 include "COMMON.IOUNITS"
641 include "COMMON.VMCPAR"
642 include "COMMON.CLASSES"
643 include "COMMON.WEIGHTDER"
644 include "COMMON.SCCOR"
645 include "COMMON.SHIELD"
647 integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock
648 double precision rri,v0ij,si
649 double precision x(nvarr)
650 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
651 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
652 double precision eps_temp(ntyp,ntyp)
653 double precision ftune_eps,ftune_epsprim
654 logical in_sc_pair_list
655 external in_sc_pair_list
656 integer itor2typ(3) /10,9,20/
657 integer ind_shield /25/
661 if (lprint) write (iout,*) "x2w: nvar",nvarr
662 if (lprint) write (iout,*) "x",(x(i),i=1,nvarr)
665 if (imask(i).eq.1 .and. i.ne.ind_shield) then
668 c write (iout,*) i, wname(i)," ",ww(i)
698 write (iout,*) "Weights:"
700 write (iout,'(a,f10.5)') wname(i),ww(i)
711 if (.not. in_sc_pair_list(iti,j)) then
712 eps(iti,j)=eps_orig(iti,j)
713 & +0.5d0*(x(ii)-eps_orig(iti,iti))
714 eps(j,iti)=eps(iti,j)
728 if (mod_tor .or. mod_fourier(nloctyp).gt.0) then
729 c write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier",
730 c & mod_fourier(nloctyp)
734 do i=-ntortyp+1,ntortyp-1
735 do j=-ntortyp+1,ntortyp-1
736 if (tor_mode.eq.0) then
738 if (mask_tor(0,j,i,iblock).eq.1) then
740 weitor(0,j,i,iblock)=x(ii)
741 else if (mask_tor(0,j,i,iblock).eq.2) then
742 do k=1,nterm(j,i,iblock)
744 v1(k,j,i,iblock)=x(ii)
746 else if (mask_tor(0,j,i,iblock).gt.2) then
747 do k=1,nterm(j,i,iblock)
749 v1(k,j,i,iblock)=x(ii)
751 do k=1,nterm(j,i,iblock)
753 v2(k,j,i,iblock)=x(ii)
758 do k=1,nterm(i,j,iblock)
759 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
760 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
761 v0ij=v0ij+si*v1(k,i,j,iblock)
764 do k=1,nlor(i,j,iblock)
765 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
768 v0(-i,-j,iblock)=v0ij
770 else if (tor_mode.eq.1) then
771 if (mask_tor(0,j,i,1).eq.1) then
773 weitor(0,j,i,1)=x(ii)
774 else if (mask_tor(0,j,i,1).eq.2) then
775 do k=1,nterm_kcc(j,i)
776 do l=1,nterm_kcc_Tb(j,i)
777 do ll=1,nterm_kcc_Tb(j,i)
779 v1_kcc(ll,l,k,j,i)=x(ii)
783 else if (mask_tor(0,j,i,1).gt.2) then
784 do k=1,nterm_kcc(j,i)
785 do l=1,nterm_kcc_Tb(j,i)
786 do ll=1,nterm_kcc_Tb(j,i)
788 v1_kcc(ll,l,k,j,i)=x(ii)
792 do k=1,nterm_kcc(j,i)
793 do l=1,nterm_kcc_Tb(j,i)
794 do ll=1,nterm_kcc_Tb(j,i)
796 v2_kcc(ll,l,k,j,i)=x(ii)
803 c Compute torsional coefficients from local energy surface expansion coeffs
804 if (mask_tor(0,j,i,1).eq.1) then
806 weitor(0,j,i,1)=x(ii)
807 c else if (mask_tor(0,j,i,1).eq.2) then
810 write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
820 do i=-nsccortyp,nsccortyp
821 do j=-nsccortyp,nsccortyp
823 if (mask_tor(k,j,i,1).eq.1) then
825 weitor(k,j,i,1)=x(ii)
826 else if (mask_tor(k,j,i,1).eq.2) then
829 v1sccor(l,k,j,i)=x(ii)
831 else if (mask_tor(k,j,i,1).gt.2) then
832 do l=1,nterm_sccor(j,i)
834 v1sccor(l,k,j,i)=x(ii)
836 do l=1,nterm_sccor(j,i)
838 v2sccor(l,k,j,i)=x(ii)
847 c 7/8/17 AL Optimization of the bending potentials
851 if (mask_ang(i).eq.0) cycle
852 do j=1,nbend_kcc_TB(i)
854 v1bend_chyb(j,i)=x(ii)
855 v1bend_chyb(j,-i)=x(ii)
861 if (mod_fourier(nloctyp).gt.0) then
866 if (mod_fourier(i).gt.0) then
870 if (mask_bnew1(k,j,i).gt.0) then
879 if (mask_bnew2(k,j,i).gt.0) then
888 if (mask_ccnew(k,j,i).gt.0) then
897 if (mask_ddnew(k,j,i).gt.0) then
905 if (mask_e0new(j,i).gt.0) then
914 if (mask_eenew(l,k,j,i).gt.0) then
928 bnew1(j,1,-i)= bnew1(j,1,i)
929 bnew1(j,2,-i)=-bnew1(j,2,i)
930 bnew2(j,1,-i)= bnew2(j,1,i)
931 bnew2(j,2,-i)=-bnew2(j,2,i)
934 c ccnew(j,1,i)=ccnew(j,1,i)/2
935 c ccnew(j,2,i)=ccnew(j,2,i)/2
936 c ddnew(j,1,i)=ddnew(j,1,i)/2
937 c ddnew(j,2,i)=ddnew(j,2,i)/2
938 ccnew(j,1,-i)=ccnew(j,1,i)
939 ccnew(j,2,-i)=-ccnew(j,2,i)
940 ddnew(j,1,-i)=ddnew(j,1,i)
941 ddnew(j,2,-i)=-ddnew(j,2,i)
943 e0new(1,-i)= e0new(1,i)
944 e0new(2,-i)=-e0new(2,i)
945 e0new(3,-i)=-e0new(3,i)
947 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
948 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
949 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
950 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
954 write (iout,'(a)') "Coefficients of the multibody terms"
955 do i=-nloctyp+1,nloctyp-1
956 write (iout,*) "Type: ",onelet(iloctyp(i))
957 write (iout,*) "Coefficients of the expansion of B1"
959 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
961 write (iout,*) "Coefficients of the expansion of B2"
963 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
965 write (iout,*) "Coefficients of the expansion of C"
966 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
967 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
968 write (iout,*) "Coefficients of the expansion of D"
969 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
970 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
971 write (iout,*) "Coefficients of the expansion of E"
972 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
975 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
982 IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN
985 if (mod_fouriertor(i).gt.0) then
989 if (mask_bnew1tor(k,j,i).gt.0) then
991 bnew1tor(k,j,i)=x(ii)
998 if (mask_bnew2tor(k,j,i).gt.0) then
1000 bnew2tor(k,j,i)=x(ii)
1007 if (mask_ccnewtor(k,j,i).gt.0) then
1009 ccnewtor(k,j,i)=x(ii)
1016 if (mask_ddnewtor(k,j,i).gt.0) then
1018 ddnewtor(k,j,i)=x(ii)
1024 if (mask_e0newtor(j,i).gt.0) then
1033 if (mask_eenewtor(l,k,j,i).gt.0) then
1035 eenewtor(l,k,j,i)=x(ii)
1048 bnew1tor(j,1,-i)= bnew1tor(j,1,i)
1049 bnew1tor(j,2,-i)=-bnew1tor(j,2,i)
1050 bnew2tor(j,1,-i)= bnew2tor(j,1,i)
1051 bnew2tor(j,2,-i)=-bnew2tor(j,2,i)
1054 c ccnew(j,1,i)=ccnew(j,1,i)/2
1055 c ccnew(j,2,i)=ccnew(j,2,i)/2
1056 c ddnew(j,1,i)=ddnew(j,1,i)/2
1057 c ddnew(j,2,i)=ddnew(j,2,i)/2
1058 ccnewtor(j,1,-i)=ccnewtor(j,1,i)
1059 ccnewtor(j,2,-i)=-ccnewtor(j,2,i)
1060 ddnewtor(j,1,-i)=ddnewtor(j,1,i)
1061 ddnewtor(j,2,-i)=-ddnewtor(j,2,i)
1063 e0newtor(1,-i)= e0newtor(1,i)
1064 e0newtor(2,-i)=-e0newtor(2,i)
1065 e0newtor(3,-i)=-e0newtor(3,i)
1067 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1068 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1069 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1070 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1075 & "Coefficients of the single-body torsional terms"
1076 do i=-nloctyp+1,nloctyp-1
1077 write (iout,*) "Type: ",onelet(iloctyp(i))
1078 write (iout,*) "Coefficients of the expansion of B1"
1080 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1082 write (iout,*) "Coefficients of the expansion of B2"
1084 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1086 write (iout,*) "Coefficients of the expansion of C"
1087 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1088 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1089 write (iout,*) "Coefficients of the expansion of D"
1090 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1091 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1092 write (iout,*) "Coefficients of the expansion of E"
1093 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1096 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1105 do i=-nloctyp+1,nloctyp-1
1108 bnew1tor(k,j,i)=bnew1(k,j,i)
1113 bnew2tor(k,j,i)=bnew2(k,j,i)
1117 ccnewtor(k,1,i)=ccnew(k,1,i)
1118 ccnewtor(k,2,i)=ccnew(k,2,i)
1119 ddnewtor(k,1,i)=ddnew(k,1,i)
1120 ddnewtor(k,2,i)=ddnew(k,2,i)
1124 ENDIF ! SPLIT_FOURIERTOR
1126 if (tor_mode.eq.2) then
1127 c Recalculate new torsional potential coefficients
1128 do i=-ntortyp+1,ntortyp-1
1129 do j=-ntortyp+1,ntortyp-1
1130 do k=1,nterm_kcc_Tb(j,i)
1131 do l=1,nterm_kcc_Tb(j,i)
1132 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1133 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1134 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1135 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1138 do k=1,nterm_kcc_Tb(j,i)
1139 do l=1,nterm_kcc_Tb(j,i)
1140 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1141 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1142 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1143 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1152 if (mod_fourier(i).gt.0) then
1155 if (mask_fourier(j,i).gt.0) then
1162 write (iout,*) "Fourier coefficients of cumulants"
1163 write (iout,*) 'Type',i
1164 write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=1,13)
1175 c B1tilde(1,i) = b(3,i)
1176 c B1tilde(2,i) =-b(5,i)
1179 CCold(1,1,i)= b(7,i)
1180 CCold(2,2,i)=-b(7,i)
1181 CCold(2,1,i)= b(9,i)
1182 CCold(1,2,i)= b(9,i)
1183 CCold(1,1,-i)= b(7,i)
1184 CCold(2,2,-i)=-b(7,i)
1185 CCold(2,1,-i)=-b(9,i)
1186 CCold(1,2,-i)=-b(9,i)
1187 c Ctilde(1,1,i)=b(7,i)
1188 c Ctilde(1,2,i)=b(9,i)
1189 c Ctilde(2,1,i)=-b(9,i)
1190 c Ctilde(2,2,i)=b(7,i)
1191 DDold(1,1,i)= b(6,i)
1192 DDold(2,2,i)=-b(6,i)
1193 DDold(2,1,i)= b(8,i)
1194 DDold(1,2,i)= b(8,i)
1195 DDold(1,1,-i)= b(6,i)
1196 DDold(2,2,-i)=-b(6,i)
1197 DDold(2,1,-i)=-b(8,i)
1198 DDold(1,2,-i)=-b(8,i)
1199 c Dtilde(1,1,i)=b(6,i)
1200 c Dtilde(1,2,i)=b(8,i)
1201 c Dtilde(2,1,i)=-b(8,i)
1202 c Dtilde(2,2,i)=b(6,i)
1203 EEold(1,1,i)= b(10,i)+b(11,i)
1204 EEold(2,2,i)=-b(10,i)+b(11,i)
1205 EEold(2,1,i)= b(12,i)-b(13,i)
1206 EEold(1,2,i)= b(12,i)+b(13,i)
1207 EEold(1,1,i)= b(10,i)+b(11,i)
1208 EEold(2,2,i)=-b(10,i)+b(11,i)
1209 EEold(2,1,i)= b(12,i)-b(13,i)
1210 EEold(1,2,i)= b(12,i)+b(13,i)
1211 EEold(1,1,-i)= b(10,i)+b(11,i)
1212 EEold(2,2,-i)=-b(10,i)+b(11,i)
1213 EEold(2,1,-i)=-b(12,i)+b(13,i)
1214 EEold(1,2,-i)=-b(12,i)-b(13,i)
1222 &"Coefficients of the cumulants (independent of valence angles)"
1223 do i=-nloctyp+1,nloctyp-1
1224 write (iout,*) 'Type ',onelet(iloctyp(i))
1226 write(iout,'(2f10.5)') B(3,i),B(5,i)
1228 write(iout,'(2f10.5)') B(2,i),B(4,i)
1231 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1235 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1239 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1246 c SC sigmas and anisotropies as parameters
1247 if (mod_side_other) then
1249 c write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4)
1252 if (mask_sigma(1,i).gt.0) then
1255 c write(iout,*) "sig0",i,ii,x(ii)
1257 if (mask_sigma(2,i).gt.0) then
1260 c write(iout,*) "sigi",i,ii,x(ii)
1262 if (mask_sigma(3,i).gt.0) then
1265 c write(iout,*) "chip",i,ii,x(ii)
1267 if (mask_sigma(4,i).gt.0) then
1270 c write(iout,*) "alp",i,ii,x(ii)
1275 if (mod_side.or.mod_side_other) then
1276 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
1277 & potname(ipot),', exponents are ',expon,2*expon
1278 C-----------------------------------------------------------------------
1279 C Calculate the "working" parameters of SC interactions.
1280 dwa16=2.0d0**(1.0d0/6.0d0)
1283 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1288 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1289 sigma(j,i)=sigma(i,j)
1290 rs0(i,j)=dwa16*sigma(i,j)
1294 if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1295 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1296 & 'Working parameters of the SC interactions:',
1297 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1302 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1311 c sigeps=dsign(1.0D0,epsij)
1312 aa(i,j)=ftune_eps(epsij)*rrij*rrij
1317 sigt1sq=sigma0(i)**2
1318 sigt2sq=sigma0(j)**2
1321 ratsig1=sigt2sq/sigt1sq
1322 ratsig2=1.0D0/ratsig1
1323 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1324 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1325 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1329 sigmaii(i,j)=rsum_max
1330 sigmaii(j,i)=rsum_max
1331 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
1333 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1334 augm(i,j)=epsij*r_augm**(2*expon)
1335 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1342 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1343 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1344 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1355 if (mask_elec(i,j,1).gt.0) then
1358 if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0
1361 if (mask_elec(i,j,2).gt.0) then
1363 rpp(i,j)=dabs(x(ii))
1364 rpp(j,i)=dabs(x(ii))
1366 if (mask_elec(i,j,3).gt.0) then
1369 if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0
1370 elpp6(j,i)=elpp6(i,j)
1372 if (mask_elec(i,j,4).gt.0) then
1378 app (i,j)=epp(i,j)*rri*rri
1379 app (j,i)=epp(j,i)*rri*rri
1380 bpp (i,j)=-2.0D0*epp(i,j)*rri
1381 bpp (j,i)=-2.0D0*epp(i,j)*rri
1382 ael6(i,j)=elpp6(i,j)*4.2D0**6
1383 ael6(j,i)=elpp6(i,j)*4.2D0**6
1384 ael3(i,j)=elpp3(i,j)*4.2D0**3
1385 ael3(j,i)=elpp3(i,j)*4.2D0**3
1389 write (iout,'(/a)') 'Electrostatic interaction constants:'
1390 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1391 & 'IT','JT','APP','BPP','AEL6','AEL3'
1394 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1395 & ael6(i,j),ael3(i,j)
1402 C Define the SC-p interaction constants
1406 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1407 & mask_scp(0,2,2) .gt. 0) then
1409 if (mask_scp(0,j,1).gt.0) then
1415 if (mask_scp(0,j,2).gt.0) then
1425 if (mask_scp(i,j,1).gt.0) then
1429 if (mask_scp(i,j,2).gt.0) then
1437 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1438 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1439 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1440 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1444 write (iout,*) "Parameters of SC-p interactions:"
1446 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1447 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1453 c Weight of the shielding term
1454 if (imask(ind_shield).eq.1) then
1456 ww(ind_shield)=x(ii)
1458 wshield=ww(ind_shield)
1460 c write (iout,*) "x2w: nvar",nvarr
1462 c Base of heat capacity
1465 if (nbeta(2,iprot).gt.0) then
1467 heatbase(iprot)=x(ii)
1471 if (ii.ne.nvarr) then
1472 write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)"
1473 write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii
1474 stop "CHUJ ABSOLUTNY!!!!"
1479 c---------------------------------------------------------------------------
1480 double precision function ftune_eps(x)
1482 double precision x,y
1483 double precision delta /1.0d-5/
1484 if (dabs(x).lt.delta) then
1486 ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8
1492 c---------------------------------------------------------------------------
1493 double precision function ftune_epsprim(x)
1495 double precision x,y
1496 double precision delta /1.0d-5/
1497 if (dabs(x).lt.delta) then
1499 ftune_epsprim=1.5d0*y-0.5*y**3
1501 ftune_epsprim=dsign(1.0d0,x)
1505 c---------------------------------------------------------------------------
1506 double precision function ftune_epsbis(x)
1508 double precision x,y
1509 double precision delta /1.0d-5/
1510 if (dabs(x).lt.delta) then
1512 ftune_epsbis=(1.5d0-1.5d0*y**3)/delta
1518 c---------------------------------------------------------------------------
1519 logical function in_sc_pair_list(iti,itj)
1521 include "DIMENSIONS"
1522 include "DIMENSIONS.ZSCOPT"
1523 include "COMMON.WEIGHTS"
1526 if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or.
1527 & ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then
1528 in_sc_pair_list=.true.
1532 in_sc_pair_list=.false.