1 subroutine w2x(nvarr,x,*)
4 include "DIMENSIONS.ZSCOPT"
7 integer IERROR,ERRCODE,kolor,key,comm
10 include "COMMON.CONTROL"
11 include "COMMON.INTERACT"
12 include "COMMON.TORSION"
13 include "COMMON.WEIGHTS"
14 include "COMMON.XBOUND"
15 include "COMMON.ENERGIES"
16 include "COMMON.IOUNITS"
17 include "COMMON.NAMES"
19 include "COMMON.VMCPAR"
20 include "COMMON.CLASSES"
21 include "COMMON.WEIGHTDER"
22 include "COMMON.FFIELD"
23 include "COMMON.SCCOR"
24 integer itor2typ(-2:2) /-20,-20,10,9,20/
26 double precision x(max_paropt),xchuj,epsi
27 integer nvarr,i,ii,j,k,l,ll,iprot,iblock
32 if (imask(i).eq.1) then
41 c Change 9/21/07 AL Variables are epsilons and not the increments on epsilons
42 c to enable control on their deviations from original values.
43 c Single-body epsilons are now defined as epsilons of pairs of same type side
46 if (me.eq.Master) then
47 write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc,
48 & " npair_sc",npair_sc
54 eps_orig(i,j)=eps(i,j)
62 c write (iout,*) "ii",ii," x",x(ii)
65 c if (epss_up(i).eq.0.0d0) then
66 c x_up(ii)=dabs(epsi)*epss_low(i)
67 c x_low(ii)=epsi-x_up(ii)
68 c x_up(ii)=epsi+x_up(ii)
70 c x_low(ii)=epss_low(i)
80 if (ityp_ssc(j).eq.iti) then
82 & write (iout,*) "Error - pair ",restyp(iti),"-",
84 & " specified in variable epsilon list but type defined",
85 & " in single-type list."
94 c if (epsp_up(i).eq.0.0d0) then
95 c x_up(ii)=dabs(epsi)*epsp_low(i)
96 c x_low(ii)=epsi-x_up(ii)
97 c x_up(ii)=epsi+x_up(ii)
99 c x_low(ii)=epsp_low(i)
100 c x_up(ii)=epsp_up(i)
102 c write (iout,*) "ii",ii," x",x(ii)
113 & write (iout,*) "Indices of optimizable torsionals"
114 do i=-ntortyp+1,ntortyp-1
115 do j=-ntortyp+1,ntortyp-1
116 c write (iout,*) "mask",i,j,mask_tor(j,i)
118 if (tor_mode.eq.0) then
121 if (mask_tor(0,j,i,iblock).eq.1) then
124 x(ii)=weitor(0,j,i,iblock)
125 x_low(ii)=weitor_low(0,j,i,iblock)
126 x_up(ii)=weitor_up(0,j,i,iblock)
128 & write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
129 & restyp(itor2typ(i)),iblock,ii
130 else if (mask_tor(0,j,i,iblock).eq.2) then
131 do k=1,nterm(j,i,iblock)
134 x(ii)=v1(k,j,i,iblock)
135 x_low(ii)=weitor_low(0,j,i,iblock)
136 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 else if (mask_tor(0,j,i,iblock).gt.2) then
142 do k=1,nterm(j,i,iblock)
145 x(ii)=v1(k,j,i,iblock)
146 x_low(ii)=weitor_low(0,j,i,iblock)
147 x_up(ii)=weitor_up(0,j,i,iblock)
149 & write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
150 & restyp(itor2typ(i)),iblock,k,ii
152 do k=1,nterm(j,i,iblock)
155 x(ii)=v2(k,j,i,iblock)
156 x_low(ii)=weitor_low(0,j,i,iblock)
157 x_up(ii)=weitor_up(0,j,i,iblock)
159 & write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
160 & restyp(itor2typ(i)),ii,k,ii
165 else if (tor_mode.eq.1) then
167 if (mask_tor(0,j,i,1).eq.1) then
170 x(ii)=weitor(0,j,i,1)
171 x_low(ii)=weitor_low(0,j,i,1)
172 x_up(ii)=weitor_up(0,j,i,1)
174 & write (iout,'(2a4,i4)') restyp(itor2typ(j)),
175 & restyp(itor2typ(i)),ii
176 else if (mask_tor(0,j,i,1).eq.2) then
177 do k=1,nterm_kcc(j,i)
178 do l=1,nterm_kcc_Tb(j,i)
179 do ll=1,nterm_kcc_Tb(j,i)
182 x(ii)=v1_kcc(ll,l,k,j,i)
183 x_low(ii)=weitor_low(0,j,i,1)
184 x_up(ii)=weitor_up(0,j,i,1)
186 & write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
187 & restyp(itor2typ(i)),k,ii
191 else if (mask_tor(0,j,i,1).gt.2) then
192 do k=1,nterm_kcc(j,i)
193 do l=1,nterm_kcc_Tb(j,i)
194 do ll=1,nterm_kcc_Tb(j,i)
197 x(ii)=v1_kcc(ll,l,k,j,i)
198 x_low(ii)=weitor_low(0,j,i,1)
199 x_up(ii)=weitor_up(0,j,i,1)
201 & write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
202 & restyp(itor2typ(i)),k,ii
206 do k=1,nterm_kcc(j,i)
207 do l=1,nterm_kcc_Tb(j,i)
208 do ll=1,nterm_kcc_Tb(j,i)
211 x(ii)=v2_kcc(ll,l,k,j,i)
212 x_low(ii)=weitor_low(0,j,i,1)
213 x_up(ii)=weitor_up(0,j,i,1)
215 & write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
216 & restyp(itor2typ(i)),k,ii
224 & write (iout,*) "ERROR! TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
232 nbacktor_var=ntor_var
233 if (me.eq.Master) write (iout,*) "nbacktor_var",nbacktor_var
238 & write (iout,*) "Indices of optimizable SCCOR torsionals"
239 do i=-nsccortyp,nsccortyp
240 do j=-nsccortyp,nsccortyp
242 c write (iout,*) "mask",i,j,mask_tor(j,i)
243 if (mask_tor(k,j,i,1).eq.1) then
246 x(ii)=weitor(k,j,i,1)
247 x_low(ii)=weitor_low(k,j,i,1)
248 x_up(ii)=weitor_up(k,j,i,1)
250 & write (iout,'(i4,2a4,i4)') k,restyp(j),restyp(i),ii
251 else if (mask_tor(k,j,i,1).eq.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)
259 & write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
261 else if (mask_tor(k,j,i,1).gt.2) then
262 do l=1,nterm_sccor(j,i)
265 x(ii)=v1sccor(l,k,j,i)
266 x_low(ii)=weitor_low(k,j,i,1)
267 x_up(ii)=weitor_up(k,j,i,1)
269 & write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
271 do l=1,nterm_sccor(j,i)
274 x(ii)=v2sccor(l,k,j,i)
275 x_low(ii)=weitor_low(k,j,i,1)
276 x_up(ii)=weitor_up(k,j,i,1)
278 & write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
287 c Optimized cumulant coefficients
289 if (mod_fourier(nloctyp).gt.0) then
294 if (mod_fourier(i).gt.0) then
298 if (mask_bnew1(k,j,i).gt.0) then
301 x_low(ii)=bnew1_low(k,j,i)
302 x_up(ii)=bnew1_up(k,j,i)
309 if (mask_bnew2(k,j,i).gt.0) then
312 x_low(ii)=bnew2_low(k,j,i)
313 x_up(ii)=bnew2_up(k,j,i)
320 if (mask_ccnew(k,j,i).gt.0) then
323 x_low(ii)=ccnew_low(k,j,i)
324 x_up(ii)=ccnew_up(k,j,i)
331 if (mask_ddnew(k,j,i).gt.0) then
334 x_low(ii)=ddnew_low(k,j,i)
335 x_up(ii)=ddnew_up(k,j,i)
341 if (mask_e0new(j,i).gt.0) then
344 x_low(ii)=e0new_low(j,i)
345 x_up(ii)=e0new_up(j,i)
352 if (mask_eenew(l,k,j,i).gt.0) then
355 x_low(ii)=eenew_low(l,k,j,i)
356 x_up(ii)=eenew_up(l,k,j,i)
368 if (mod_fourier(i).gt.0) then
370 if (mask_fourier(j,i).gt.0) then
383 c Added SC sigmas and anisotropies
385 if (mod_side_other) then
387 if (mask_sigma(1,i).gt.0) then
390 x_low(ii)=sigma_low(1,i)
391 x_up(ii)=sigma_up(1,i)
392 c write (iout,*) "sig0",i,ii,x(ii), x_low(ii),x_up(ii)
394 if (mask_sigma(2,i).gt.0) then
397 x_low(ii)=sigma_low(2,i)
398 x_up(ii)=sigma_up(2,i)
399 c write (iout,*) "sigi",i,ii,x(ii), x_low(ii),x_up(ii)
401 if (mask_sigma(3,i).gt.0) then
404 x_low(ii)=sigma_low(3,i)
405 x_up(ii)=sigma_up(3,i)
406 c write (iout,*) "chi",i,ii,x(ii), x_low(ii),x_up(ii)
408 if (mask_sigma(4,i).gt.0) then
411 x_low(ii)=sigma_low(4,i)
412 x_up(ii)=sigma_up(4,i)
413 c write (iout,*) "alp",i,ii,x(ii), x_low(ii),x_up(ii)
423 if (mask_elec(i,j,1).gt.0) then
426 x_low(ii)=epp_low(i,j)
429 if (mask_elec(i,j,2).gt.0) then
432 x_low(ii)=rpp_low(i,j)
435 if (mask_elec(i,j,3).gt.0) then
438 x_low(ii)=elpp6_low(i,j)
439 x_up(ii)=elpp6_up(i,j)
441 if (mask_elec(i,j,4).gt.0) then
444 x_low(ii)=elpp3_low(i,j)
445 x_up(ii)=elpp3_up(i,j)
452 C Define the SC-p interaction constants
456 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
457 & mask_scp(0,2,2) .gt. 0) then
459 if (mask_scp(0,j,1).gt.0) then
462 x_low(ii)=epscp_low(0,j)
463 x_up(ii)=epscp_up(0,j)
465 if (mask_scp(0,j,2).gt.0) then
468 x_low(ii)=rscp_low(0,j)
469 x_up(ii)=rscp_up(0,j)
475 if (mask_scp(i,j,1).gt.0) then
478 x_low(ii)=epscp_low(i,j)
479 x_up(ii)=epscp_up(i,j)
481 if (mask_scp(i,j,2).gt.0) then
484 x_low(ii)=rscp_low(i,j)
485 x_up(ii)=rscp_up(i,j)
494 c Compute boundaries, if defined relative to starting values
497 & write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i)
498 if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then
502 else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then
504 if (x(i).gt.0.0d0) then
505 x_low(i)=x(i)*(1.0d0-xchuj)
506 x_up(i)=x(i)*(1.0d0+xchuj)
508 x_low(i)=x(i)*(1.0d0+xchuj)
509 x_up(i)=x(i)*(1.0d0-xchuj)
514 c Base of heat capacity
517 if (nbeta(2,iprot).gt.0) then
519 x(ii)=heatbase(iprot)
527 if (me.eq.Master)then
528 write (iout,*) "Number of optimizable parameters:",nvarr
530 write (iout,*) "Initial variables and boundaries:"
531 write (iout,'(3x,3a10)') "x","x_low","x_up"
533 write (iout,'(i3,3f10.5)') i,x(i),x_low(i),x_up(i)
539 c------------------------------------------------------------------------------
540 subroutine x2w(nvarr,x)
543 include "DIMENSIONS.ZSCOPT"
546 integer IERROR,ERRCODE,kolor,key,comm
549 include "COMMON.CONTROL"
550 include "COMMON.NAMES"
551 include "COMMON.WEIGHTS"
552 include "COMMON.XBOUND"
553 include "COMMON.INTERACT"
554 include "COMMON.ENERGIES"
555 include "COMMON.FFIELD"
556 include "COMMON.TORSION"
557 include "COMMON.IOUNITS"
558 include "COMMON.VMCPAR"
559 include "COMMON.CLASSES"
560 include "COMMON.WEIGHTDER"
561 include "COMMON.SCCOR"
563 integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock
564 double precision rri,v0ij,si
565 double precision x(nvarr)
566 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
567 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
568 double precision eps_temp(ntyp,ntyp)
569 double precision ftune_eps,ftune_epsprim
570 logical in_sc_pair_list
571 external in_sc_pair_list
572 integer itor2typ(3) /10,9,20/
575 if (me.ne.Master) lprint = .false.
577 if (lprint) write (iout,*) "x2w: nvar",nvarr
578 if (lprint) write (iout,*) "x",(x(i),i=1,nvarr)
581 if (imask(i).eq.1) then
613 write (iout,*) "Weights:"
615 write (iout,'(a,f10.5)') wname(i),ww(i)
626 if (.not. in_sc_pair_list(iti,j)) then
627 eps(iti,j)=eps_orig(iti,j)
628 & +0.5d0*(x(ii)-eps_orig(iti,iti))
629 eps(j,iti)=eps(iti,j)
643 if (mod_tor .or. mod_fourier(nloctyp).gt.0) then
644 c write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier",
645 c & mod_fourier(nloctyp)
649 do i=-ntortyp+1,ntortyp-1
650 do j=-ntortyp+1,ntortyp-1
651 if (tor_mode.eq.0) then
653 if (mask_tor(0,j,i,iblock).eq.1) then
655 weitor(0,j,i,iblock)=x(ii)
656 else if (mask_tor(0,j,i,iblock).eq.2) then
657 do k=1,nterm(j,i,iblock)
659 v1(k,j,i,iblock)=x(ii)
661 else if (mask_tor(0,j,i,iblock).gt.2) then
662 do k=1,nterm(j,i,iblock)
664 v1(k,j,i,iblock)=x(ii)
666 do k=1,nterm(j,i,iblock)
668 v2(k,j,i,iblock)=x(ii)
673 do k=1,nterm(i,j,iblock)
674 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
675 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
676 v0ij=v0ij+si*v1(k,i,j,iblock)
679 do k=1,nlor(i,j,iblock)
680 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
683 v0(-i,-j,iblock)=v0ij
685 else if (tor_mode.eq.1) then
686 if (mask_tor(0,j,i,1).eq.1) then
688 weitor(0,j,i,1)=x(ii)
689 else if (mask_tor(0,j,i,1).eq.2) then
690 do k=1,nterm_kcc(j,i)
691 do l=1,nterm_kcc_Tb(j,i)
692 do ll=1,nterm_kcc_Tb(j,i)
694 v1_kcc(ll,l,k,j,i)=x(ii)
698 else if (mask_tor(0,j,i,1).gt.2) then
699 do k=1,nterm_kcc(j,i)
700 do l=1,nterm_kcc_Tb(j,i)
701 do ll=1,nterm_kcc_Tb(j,i)
703 v1_kcc(ll,l,k,j,i)=x(ii)
707 do k=1,nterm_kcc(j,i)
708 do l=1,nterm_kcc_Tb(j,i)
709 do ll=1,nterm_kcc_Tb(j,i)
711 v2_kcc(ll,l,k,j,i)=x(ii)
718 c Compute torsional coefficients from local energy surface expansion coeffs
719 if (mask_tor(0,j,i,1).eq.1) then
721 weitor(0,j,i,1)=x(ii)
722 c else if (mask_tor(0,j,i,1).eq.2) then
726 write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
736 do i=-nsccortyp,nsccortyp
737 do j=-nsccortyp,nsccortyp
739 if (mask_tor(k,j,i,1).eq.1) then
741 weitor(k,j,i,1)=x(ii)
742 else if (mask_tor(k,j,i,1).eq.2) then
745 v1sccor(l,k,j,i)=x(ii)
747 else if (mask_tor(k,j,i,1).gt.2) then
748 do l=1,nterm_sccor(j,i)
750 v1sccor(l,k,j,i)=x(ii)
752 do l=1,nterm_sccor(j,i)
754 v2sccor(l,k,j,i)=x(ii)
763 if (mod_fourier(nloctyp).gt.0) then
768 if (mod_fourier(i).gt.0) then
772 if (mask_bnew1(k,j,i).gt.0) then
781 if (mask_bnew2(k,j,i).gt.0) then
790 if (mask_ccnew(k,j,i).gt.0) then
799 if (mask_ddnew(k,j,i).gt.0) then
807 if (mask_e0new(j,i).gt.0) then
816 if (mask_eenew(l,k,j,i).gt.0) then
830 bnew1(j,1,-i)= bnew1(j,1,i)
831 bnew1(j,2,-i)=-bnew1(j,2,i)
832 bnew2(j,1,-i)= bnew2(j,1,i)
833 bnew2(j,2,-i)=-bnew2(j,2,i)
836 c ccnew(j,1,i)=ccnew(j,1,i)/2
837 c ccnew(j,2,i)=ccnew(j,2,i)/2
838 c ddnew(j,1,i)=ddnew(j,1,i)/2
839 c ddnew(j,2,i)=ddnew(j,2,i)/2
840 ccnew(j,1,-i)=ccnew(j,1,i)
841 ccnew(j,2,-i)=-ccnew(j,2,i)
842 ddnew(j,1,-i)=ddnew(j,1,i)
843 ddnew(j,2,-i)=-ddnew(j,2,i)
845 e0new(1,-i)= e0new(1,i)
846 e0new(2,-i)=-e0new(2,i)
847 e0new(3,-i)=-e0new(3,i)
849 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
850 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
851 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
852 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
856 write (iout,'(a)') "Coefficients of the multibody terms"
857 do i=-nloctyp+1,nloctyp-1
858 write (iout,*) "Type: ",onelet(iloctyp(i))
859 write (iout,*) "Coefficients of the expansion of B1"
861 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
863 write (iout,*) "Coefficients of the expansion of B2"
865 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
867 write (iout,*) "Coefficients of the expansion of C"
868 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
869 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
870 write (iout,*) "Coefficients of the expansion of D"
871 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
872 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
873 write (iout,*) "Coefficients of the expansion of E"
874 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
877 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
883 if (tor_mode.eq.2) then
884 c Recalculate new torsional potential coefficients
885 do i=-ntortyp+1,ntortyp-1
886 do j=-ntortyp+1,ntortyp-1
887 do k=1,nterm_kcc_Tb(j,i)
888 do l=1,nterm_kcc_Tb(j,i)
889 v1_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,1,j)
890 & +bnew1(k,2,i)*bnew2(l,2,j)
891 v2_kcc(k,l,1,i,j)=bnew1(k,1,i)*bnew2(l,2,j)
892 & +bnew1(k,2,i)*bnew2(l,1,j)
895 do k=1,nterm_kcc_Tb(j,i)
896 do l=1,nterm_kcc_Tb(j,i)
897 v1_kcc(k,l,2,i,j)=-0.25*(ccnew(k,1,i)*ddnew(l,1,j)
898 & -ccnew(k,2,i)*ddnew(l,2,j))
899 v2_kcc(k,l,2,i,j)=-0.25*(ccnew(k,2,i)*ddnew(l,1,j)
900 & +ccnew(k,1,i)*ddnew(l,2,j))
909 if (mod_fourier(i).gt.0) then
912 if (mask_fourier(j,i).gt.0) then
919 write (iout,*) "Fourier coefficients of cumulants"
920 write (iout,*) 'Type',i
921 write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=1,13)
932 c B1tilde(1,i) = b(3,i)
933 c B1tilde(2,i) =-b(5,i)
940 CCold(1,1,-i)= b(7,i)
941 CCold(2,2,-i)=-b(7,i)
942 CCold(2,1,-i)=-b(9,i)
943 CCold(1,2,-i)=-b(9,i)
944 c Ctilde(1,1,i)=b(7,i)
945 c Ctilde(1,2,i)=b(9,i)
946 c Ctilde(2,1,i)=-b(9,i)
947 c Ctilde(2,2,i)=b(7,i)
952 DDold(1,1,-i)= b(6,i)
953 DDold(2,2,-i)=-b(6,i)
954 DDold(2,1,-i)=-b(8,i)
955 DDold(1,2,-i)=-b(8,i)
956 c Dtilde(1,1,i)=b(6,i)
957 c Dtilde(1,2,i)=b(8,i)
958 c Dtilde(2,1,i)=-b(8,i)
959 c Dtilde(2,2,i)=b(6,i)
960 EEold(1,1,i)= b(10,i)+b(11,i)
961 EEold(2,2,i)=-b(10,i)+b(11,i)
962 EEold(2,1,i)= b(12,i)-b(13,i)
963 EEold(1,2,i)= b(12,i)+b(13,i)
964 EEold(1,1,i)= b(10,i)+b(11,i)
965 EEold(2,2,i)=-b(10,i)+b(11,i)
966 EEold(2,1,i)= b(12,i)-b(13,i)
967 EEold(1,2,i)= b(12,i)+b(13,i)
968 EEold(1,1,-i)= b(10,i)+b(11,i)
969 EEold(2,2,-i)=-b(10,i)+b(11,i)
970 EEold(2,1,-i)=-b(12,i)+b(13,i)
971 EEold(1,2,-i)=-b(12,i)-b(13,i)
979 &"Coefficients of the cumulants (independent of valence angles)"
980 do i=-nloctyp+1,nloctyp-1
981 write (iout,*) 'Type ',onelet(iloctyp(i))
983 write(iout,'(2f10.5)') B(3,i),B(5,i)
985 write(iout,'(2f10.5)') B(2,i),B(4,i)
988 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
992 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
996 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1003 c SC sigmas and anisotropies as parameters
1004 if (mod_side_other) then
1006 c write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4)
1009 if (mask_sigma(1,i).gt.0) then
1012 c write(iout,*) "sig0",i,ii,x(ii)
1014 if (mask_sigma(2,i).gt.0) then
1017 c write(iout,*) "sigi",i,ii,x(ii)
1019 if (mask_sigma(3,i).gt.0) then
1022 c write(iout,*) "chip",i,ii,x(ii)
1024 if (mask_sigma(4,i).gt.0) then
1027 c write(iout,*) "alp",i,ii,x(ii)
1032 if (mod_side.or.mod_side_other) then
1033 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
1034 & potname(ipot),', exponents are ',expon,2*expon
1035 C-----------------------------------------------------------------------
1036 C Calculate the "working" parameters of SC interactions.
1037 dwa16=2.0d0**(1.0d0/6.0d0)
1040 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1045 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1046 sigma(j,i)=sigma(i,j)
1047 rs0(i,j)=dwa16*sigma(i,j)
1051 if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1052 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1053 & 'Working parameters of the SC interactions:',
1054 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1059 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1068 c sigeps=dsign(1.0D0,epsij)
1069 aa(i,j)=ftune_eps(epsij)*rrij*rrij
1074 sigt1sq=sigma0(i)**2
1075 sigt2sq=sigma0(j)**2
1078 ratsig1=sigt2sq/sigt1sq
1079 ratsig2=1.0D0/ratsig1
1080 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1081 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1082 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1086 sigmaii(i,j)=rsum_max
1087 sigmaii(j,i)=rsum_max
1088 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
1090 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1091 augm(i,j)=epsij*r_augm**(2*expon)
1092 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1099 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1100 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1101 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1112 if (mask_elec(i,j,1).gt.0) then
1115 if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0
1118 if (mask_elec(i,j,2).gt.0) then
1120 rpp(i,j)=dabs(x(ii))
1121 rpp(j,i)=dabs(x(ii))
1123 if (mask_elec(i,j,3).gt.0) then
1126 if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0
1127 elpp6(j,i)=elpp6(i,j)
1129 if (mask_elec(i,j,4).gt.0) then
1135 app (i,j)=epp(i,j)*rri*rri
1136 app (j,i)=epp(j,i)*rri*rri
1137 bpp (i,j)=-2.0D0*epp(i,j)*rri
1138 bpp (j,i)=-2.0D0*epp(i,j)*rri
1139 ael6(i,j)=elpp6(i,j)*4.2D0**6
1140 ael6(j,i)=elpp6(i,j)*4.2D0**6
1141 ael3(i,j)=elpp3(i,j)*4.2D0**3
1142 ael3(j,i)=elpp3(i,j)*4.2D0**3
1146 write (iout,'(/a)') 'Electrostatic interaction constants:'
1147 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1148 & 'IT','JT','APP','BPP','AEL6','AEL3'
1151 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1152 & ael6(i,j),ael3(i,j)
1159 C Define the SC-p interaction constants
1163 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1164 & mask_scp(0,2,2) .gt. 0) then
1166 if (mask_scp(0,j,1).gt.0) then
1172 if (mask_scp(0,j,2).gt.0) then
1182 if (mask_scp(i,j,1).gt.0) then
1186 if (mask_scp(i,j,2).gt.0) then
1194 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1195 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1196 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1197 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1201 write (iout,*) "Parameters of SC-p interactions:"
1203 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1204 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1211 c write (iout,*) "x2w: nvar",nvarr
1213 c Base of heat capacity
1216 if (nbeta(2,iprot).gt.0) then
1218 heatbase(iprot)=x(ii)
1222 if (ii.ne.nvarr) then
1223 write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)"
1224 write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii
1225 stop "CHUJ ABSOLUTNY!!!!"
1230 c---------------------------------------------------------------------------
1231 double precision function ftune_eps(x)
1233 double precision x,y
1234 double precision delta /1.0d-5/
1235 if (dabs(x).lt.delta) then
1237 ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8
1243 c---------------------------------------------------------------------------
1244 double precision function ftune_epsprim(x)
1246 double precision x,y
1247 double precision delta /1.0d-5/
1248 if (dabs(x).lt.delta) then
1250 ftune_epsprim=1.5d0*y-0.5*y**3
1252 ftune_epsprim=dsign(1.0d0,x)
1256 c---------------------------------------------------------------------------
1257 double precision function ftune_epsbis(x)
1259 double precision x,y
1260 double precision delta /1.0d-5/
1261 if (dabs(x).lt.delta) then
1263 ftune_epsbis=(1.5d0-1.5d0*y**3)/delta
1269 c---------------------------------------------------------------------------
1270 logical function in_sc_pair_list(iti,itj)
1272 include "DIMENSIONS"
1273 include "DIMENSIONS.ZSCOPT"
1274 include "COMMON.WEIGHTS"
1277 if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or.
1278 & ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then
1279 in_sc_pair_list=.true.
1283 in_sc_pair_list=.false.