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
26 c write (iout,*) "W2X: SPLIT_FOURIERTOR",SPLIT_FOURIERTOR,
27 c & " tor_mode",tor_mode
30 if (imask(i).eq.1 .and. i.ne.ind_shield) then
38 c 2/10/208 AL: insert the base torsional PMFs
39 if (mod_fourier(nloctyp).gt.0) then
45 x_low(ii)=e0_low(iti,itj)
46 x_up(ii)=e0_up(iti,itj)
57 if (torsion_pmf.and.eello_pmf) then
60 x_low(ii)=wello_PMF_low
66 if (turn_pmf.and.(torsion_pmf.or.eello_pmf)) then
69 x_low(ii)=wturn_PMF_low
80 c Change 9/21/07 AL Variables are epsilons and not the increments on epsilons
81 c to enable control on their deviations from original values.
82 c Single-body epsilons are now defined as epsilons of pairs of same type side
85 c write (iout,*) "W2X: mod_side",mod_side," nsingle_sc",nsingle_sc,
86 c & " npair_sc",npair_sc
91 eps_orig(i,j)=eps(i,j)
99 c write (iout,*) "ii",ii," x",x(ii)
100 x_low(ii)=epss_low(i)
102 c if (epss_up(i).eq.0.0d0) then
103 c x_up(ii)=dabs(epsi)*epss_low(i)
104 c x_low(ii)=epsi-x_up(ii)
105 c x_up(ii)=epsi+x_up(ii)
107 c x_low(ii)=epss_low(i)
108 c x_up(ii)=epss_up(i)
117 if (ityp_ssc(j).eq.iti) then
118 write (iout,*) "Error - pair ",restyp(iti),"-",
120 & " specified in variable epsilon list but type defined",
121 & " in single-type list."
128 x_low(ii)=epsp_low(i)
130 c if (epsp_up(i).eq.0.0d0) then
131 c x_up(ii)=dabs(epsi)*epsp_low(i)
132 c x_low(ii)=epsi-x_up(ii)
133 c x_up(ii)=epsi+x_up(ii)
135 c x_low(ii)=epsp_low(i)
136 c x_up(ii)=epsp_up(i)
138 c write (iout,*) "ii",ii," x",x(ii)
148 c write (iout,*) "Indices of optimizable torsionals"
149 do i=-ntortyp+1,ntortyp-1
150 do j=-ntortyp+1,ntortyp-1
151 c write (iout,*) "mask",i,j,mask_tor(j,i)
153 if (tor_mode.eq.0) then
156 if (mask_tor(0,j,i,iblock).eq.1) then
159 x(ii)=weitor(0,j,i,iblock)
160 x_low(ii)=weitor_low(0,j,i,iblock)
161 x_up(ii)=weitor_up(0,j,i,iblock)
162 c write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
163 c & restyp(itor2typ(i)),iblock,ii
164 else if (mask_tor(0,j,i,iblock).eq.2) then
165 do k=1,nterm(j,i,iblock)
168 x(ii)=v1(k,j,i,iblock)
169 x_low(ii)=weitor_low(0,j,i,iblock)
170 x_up(ii)=weitor_up(0,j,i,iblock)
171 c write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
172 c & restyp(itor2typ(i)),iblock,k,ii
174 else if (mask_tor(0,j,i,iblock).gt.2) then
175 do k=1,nterm(j,i,iblock)
178 x(ii)=v1(k,j,i,iblock)
179 x_low(ii)=weitor_low(0,j,i,iblock)
180 x_up(ii)=weitor_up(0,j,i,iblock)
181 c write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
182 c & restyp(itor2typ(i)),iblock,k,ii
184 do k=1,nterm(j,i,iblock)
187 x(ii)=v2(k,j,i,iblock)
188 x_low(ii)=weitor_low(0,j,i,iblock)
189 x_up(ii)=weitor_up(0,j,i,iblock)
190 c write (iout,'(2a4,3i4)') restyp(itor2typ(j)),
191 c & restyp(itor2typ(i)),ii,k,ii
196 else if (tor_mode.eq.1) then
198 if (mask_tor(0,j,i,1).eq.1) then
201 x(ii)=weitor(0,j,i,1)
202 x_low(ii)=weitor_low(0,j,i,1)
203 x_up(ii)=weitor_up(0,j,i,1)
204 c write (iout,'(2a4,i4)') restyp(itor2typ(j)),
205 c & restyp(itor2typ(i)),ii
206 else if (mask_tor(0,j,i,1).eq.2) then
207 do k=1,nterm_kcc(j,i)
208 do l=1,nterm_kcc_Tb(j,i)
209 do ll=1,nterm_kcc_Tb(j,i)
212 x(ii)=v1_kcc(ll,l,k,j,i)
213 x_low(ii)=weitor_low(0,j,i,1)
214 x_up(ii)=weitor_up(0,j,i,1)
215 c write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
216 c & restyp(itor2typ(i)),k,ii
220 else if (mask_tor(0,j,i,1).gt.2) then
221 do k=1,nterm_kcc(j,i)
222 do l=1,nterm_kcc_Tb(j,i)
223 do ll=1,nterm_kcc_Tb(j,i)
226 x(ii)=v1_kcc(ll,l,k,j,i)
227 x_low(ii)=weitor_low(0,j,i,1)
228 x_up(ii)=weitor_up(0,j,i,1)
229 c write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
230 c & restyp(itor2typ(i)),k,ii
234 do k=1,nterm_kcc(j,i)
235 do l=1,nterm_kcc_Tb(j,i)
236 do ll=1,nterm_kcc_Tb(j,i)
239 x(ii)=v2_kcc(ll,l,k,j,i)
240 x_low(ii)=weitor_low(0,j,i,1)
241 x_up(ii)=weitor_up(0,j,i,1)
242 c write (iout,'(2a4,2i4)') restyp(itor2typ(j)),
243 c & restyp(itor2typ(i)),k,ii
250 c Compute torsional coefficients from local energy surface expansion coeffs
251 if (mask_tor(0,j,i,1).eq.1) then
253 c write (iout,*) "i",i," j",j," weitor",weitor(0,j,i,1)
254 x(ii)=weitor(0,j,i,1)
255 x_low(ii)=weitor_low(0,j,i,1)
256 x_up(ii)=weitor_up(0,j,i,1)
259 write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
268 nbacktor_var=ntor_var
269 c write (iout,*) "nbacktor_var",nbacktor_var
273 c write (iout,*) "Indices of optimizable SCCOR torsionals"
274 do i=-nsccortyp,nsccortyp
275 do j=-nsccortyp,nsccortyp
277 c write (iout,*) "mask",i,j,mask_tor(j,i)
278 if (mask_tor(k,j,i,1).eq.1) then
281 x(ii)=weitor(k,j,i,1)
282 x_low(ii)=weitor_low(k,j,i,1)
283 x_up(ii)=weitor_up(k,j,i,1)
284 c write (iout,'(i4,2a4,i4)') k,restyp(j),restyp(i),ii
285 else if (mask_tor(k,j,i,1).eq.2) then
286 do l=1,nterm_sccor(j,i)
289 x(ii)=v1sccor(l,k,j,i)
290 x_low(ii)=weitor_low(k,j,i,1)
291 x_up(ii)=weitor_up(k,j,i,1)
292 c write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
294 else if (mask_tor(k,j,i,1).gt.2) then
295 do l=1,nterm_sccor(j,i)
298 x(ii)=v1sccor(l,k,j,i)
299 x_low(ii)=weitor_low(k,j,i,1)
300 x_up(ii)=weitor_up(k,j,i,1)
301 c write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
303 do l=1,nterm_sccor(j,i)
306 x(ii)=v2sccor(l,k,j,i)
307 x_low(ii)=weitor_low(k,j,i,1)
308 x_up(ii)=weitor_up(k,j,i,1)
309 c write (iout,'(i4,2a4,2i4)') k,restyp(j),restyp(i),l,ii
318 c 7/8/17 AL Optimization of the bending potentials
323 if (mask_ang(i).eq.0) cycle
324 do j=1,nbend_kcc_TB(i)
327 x(ii)=v1bend_chyb(j,i)
328 x_low(ii)=v1bend_low(j,i)
329 x_up(ii)=v1bend_up(j,i)
333 c write (iout,*) "nang_var",nang_var
335 c Optimized cumulant coefficients
337 if (mod_fourier(nloctyp).gt.0) then
342 if (mod_fourier(i).gt.0) then
346 if (mask_bnew1(k,j,i).gt.0) then
349 if (mask_bnewtor1(k,j,i).eq.0) mask_bnewtor1(k,j,i)=-ii
351 x_low(ii)=bnew1_low(k,j,i)
352 x_up(ii)=bnew1_up(k,j,i)
359 if (mask_bnew2(k,j,i).gt.0) then
363 x_low(ii)=bnew2_low(k,j,i)
364 x_up(ii)=bnew2_up(k,j,i)
371 if (mask_ccnew(k,j,i).gt.0) then
375 x_low(ii)=ccnew_low(k,j,i)
376 x_up(ii)=ccnew_up(k,j,i)
383 if (mask_ddnew(k,j,i).gt.0) then
387 x_low(ii)=ddnew_low(k,j,i)
388 x_up(ii)=ddnew_up(k,j,i)
394 if (mask_e0new(j,i).gt.0) then
398 x_low(ii)=e0new_low(j,i)
399 x_up(ii)=e0new_up(j,i)
406 if (mask_eenew(l,k,j,i).gt.0) then
408 mask_eenew(l,k,j,i)=ii
410 x_low(ii)=eenew_low(l,k,j,i)
411 x_up(ii)=eenew_up(l,k,j,i)
421 IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN
425 c write (iout,*) "type",i," mod_fouriertor",mod_fouriertor(i)
426 if (mod_fouriertor(i).gt.0) then
430 c write (iout,*) "mask_bnew2tor",k,j,mask_bnew1tor(k,j,i)
431 if (mask_bnew1tor(k,j,i).gt.0) then
433 mask_bnew1tor(k,j,i)=ii
434 x(ii)=bnew1tor(k,j,i)
435 x_low(ii)=bnew1tor_low(k,j,i)
436 x_up(ii)=bnew1tor_up(k,j,i)
443 c write (iout,*) "mask_bnew2tor",j,k,mask_bnew2tor(k,j,i)
444 if (mask_bnew2tor(k,j,i).gt.0) then
446 mask_bnew2tor(k,j,i)=ii
447 x(ii)=bnew2tor(k,j,i)
448 x_low(ii)=bnew2tor_low(k,j,i)
449 x_up(ii)=bnew2tor_up(k,j,i)
456 c write (iout,*) "mask_ccnewtor",k,j,mask_ccnewtor(k,j,i)
457 if (mask_ccnewtor(k,j,i).gt.0) then
459 mask_ccnewtor(k,j,i)=ii
460 x(ii)=ccnewtor(k,j,i)
461 x_low(ii)=ccnewtor_low(k,j,i)
462 x_up(ii)=ccnewtor_up(k,j,i)
469 c write (iout,*) "mask_ddnewtor",k,j,mask_ddnewtor(k,j,i)
470 if (mask_ddnewtor(k,j,i).gt.0) then
472 mask_ddnewtor(k,j,i)=ii
473 x(ii)=ddnewtor(k,j,i)
474 x_low(ii)=ddnewtor_low(k,j,i)
475 x_up(ii)=ddnewtor_up(k,j,i)
480 c write (iout,*) "type",i
482 c write (iout,*) "mask_e0newtor",j,mask_e0newtor(j,i)
483 if (mask_e0newtor(j,i).gt.0) then
485 mask_e0newtor(j,i)=ii
487 x_low(ii)=e0newtor_low(j,i)
488 x_up(ii)=e0newtor_up(j,i)
489 c write (iout,*) j," e0newtorlow",e0newtor_low(j,i),
490 c & " e0newtorup",e0newtor_up(j,i)
497 c write (iout,*)"mask_eenewtor",l,k,j,mask_eenewtor(l,k,j,i)
498 if (mask_eenewtor(l,k,j,i).gt.0) then
500 mask_eenewtor(l,k,j,i)=ii
501 x(ii)=eenewtor(l,k,j,i)
502 x_low(ii)=eenewtor_low(l,k,j,i)
503 x_up(ii)=eenewtor_up(l,k,j,i)
504 c write (iout,*) j,k,l,
505 c & " eenewtorlow",eenewtor_low(l,k,j,i),
506 c & " eenewtorup",eenewtor_up(l,k,j,i)
521 mask_bnewtor1(k,j,i)=-mask_bnew1(k,j,i)
527 mask_bnewtor2(k,j,i)=-mask_bnew2(k,j,i)
533 mask_ccnewtor(k,j,i)=-mask_ccnew(k,j,i)
539 mask_ddnewtor(k,j,i)=-mask_ddnew(k,j,i)
544 mask_e0newtor(j,i)=-mask_e0new(j,i)
550 mask_eenewtor(l,k,j,i)=-mask_eenew(l,k,j,i)
560 if (mod_fourier(i).gt.0) then
562 if (mask_fourier(j,i).gt.0) then
575 c Added SC sigmas and anisotropies
577 if (mod_side_other) then
579 if (mask_sigma(1,i).gt.0) then
582 x_low(ii)=sigma_low(1,i)
583 x_up(ii)=sigma_up(1,i)
584 c write (iout,*) "sig0",i,ii,x(ii), x_low(ii),x_up(ii)
586 if (mask_sigma(2,i).gt.0) then
589 x_low(ii)=sigma_low(2,i)
590 x_up(ii)=sigma_up(2,i)
591 c write (iout,*) "sigi",i,ii,x(ii), x_low(ii),x_up(ii)
593 if (mask_sigma(3,i).gt.0) then
596 x_low(ii)=sigma_low(3,i)
597 x_up(ii)=sigma_up(3,i)
598 c write (iout,*) "chi",i,ii,x(ii), x_low(ii),x_up(ii)
600 if (mask_sigma(4,i).gt.0) then
603 x_low(ii)=sigma_low(4,i)
604 x_up(ii)=sigma_up(4,i)
605 c write (iout,*) "alp",i,ii,x(ii), x_low(ii),x_up(ii)
615 if (mask_elec(i,j,1).gt.0) then
618 x_low(ii)=epp_low(i,j)
621 if (mask_elec(i,j,2).gt.0) then
624 x_low(ii)=rpp_low(i,j)
627 if (mask_elec(i,j,3).gt.0) then
630 x_low(ii)=elpp6_low(i,j)
631 x_up(ii)=elpp6_up(i,j)
633 if (mask_elec(i,j,4).gt.0) then
636 x_low(ii)=elpp3_low(i,j)
637 x_up(ii)=elpp3_up(i,j)
644 C Define the SC-p interaction constants
648 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
649 & mask_scp(0,2,2) .gt. 0) then
651 if (mask_scp(0,j,1).gt.0) then
654 x_low(ii)=epscp_low(0,j)
655 x_up(ii)=epscp_up(0,j)
657 if (mask_scp(0,j,2).gt.0) then
660 x_low(ii)=rscp_low(0,j)
661 x_up(ii)=rscp_up(0,j)
667 if (mask_scp(i,j,1).gt.0) then
670 x_low(ii)=epscp_low(i,j)
671 x_up(ii)=epscp_up(i,j)
673 if (mask_scp(i,j,2).gt.0) then
676 x_low(ii)=rscp_low(i,j)
677 x_up(ii)=rscp_up(i,j)
685 c Weight of the shielding term
686 if (imask(ind_shield).eq.1) then
689 x_low(ii)=ww_low(ind_shield)
690 x_up(ii)=ww_up(ind_shield)
694 c Compute boundaries, if defined relative to starting values
696 write (iout,*) "i",i," x_low",x_low(i)," x_up",x_up(i)
697 if (x_up(i).eq.0.0d0 .and. x_low(i).gt.0.0d0) then
701 else if (x_up(i).eq.-1.0d0 .and. x_low(i).gt.0.0d0) then
703 if (x(i).gt.0.0d0) then
704 x_low(i)=x(i)*(1.0d0-xchuj)
705 x_up(i)=x(i)*(1.0d0+xchuj)
707 x_low(i)=x(i)*(1.0d0+xchuj)
708 x_up(i)=x(i)*(1.0d0-xchuj)
713 c Base of heat capacity
716 if (nbeta(2,iprot).gt.0) then
718 x(ii)=heatbase(iprot)
726 write (iout,*) "Number of optimizable parameters:",nvarr
728 write (iout,*) "Initial variables and boundaries:"
729 write (iout,'(3x,3a15)') " x"," x_low"," x_up"
731 write (iout,'(i3,3e15.5)') i,x(i),x_low(i),x_up(i)
736 c------------------------------------------------------------------------------
737 subroutine x2w(nvarr,x)
740 include "DIMENSIONS.ZSCOPT"
741 include "COMMON.CONTROL"
742 include "COMMON.NAMES"
743 include "COMMON.WEIGHTS"
744 include "COMMON.INTERACT"
745 include "COMMON.ENERGIES"
746 include "COMMON.FFIELD"
747 include "COMMON.TORSION"
748 include "COMMON.LOCAL"
749 include "COMMON.IOUNITS"
750 include "COMMON.VMCPAR"
751 include "COMMON.CLASSES"
752 include "COMMON.WEIGHTDER"
753 include "COMMON.SCCOR"
754 include "COMMON.SHIELD"
756 integer i,j,k,kk,l,ll,ii,nvarr,iti,itj,iprot,iblock
757 double precision rri,v0ij,si
758 double precision x(nvarr)
759 double precision epsij,rrij,sigeps,sigt1sq,sigt2sq,sigii1,sigii2,
760 & ratsig1,ratsig2,rsum_max,dwa16,r_augm
761 double precision eps_temp(ntyp,ntyp)
762 double precision ftune_eps,ftune_epsprim
763 logical in_sc_pair_list
764 external in_sc_pair_list
765 integer itor2typ(3) /10,9,20/
766 integer ind_shield /25/
770 if (lprint) write (iout,*) "x2w: nvar",nvarr
771 if (lprint) write (iout,*) "x",(x(i),i=1,nvarr)
772 if (lprint) write(iout,*)"X2W: SPLIT_FOURIERTOR",SPLIT_FOURIERTOR,
773 & " tor_mode",tor_mode
776 if (imask(i).eq.1 .and. i.ne.ind_shield) then
779 c write (iout,*) i, wname(i)," ",ww(i)
783 c 2/10/208 AL: insert the base torsional PMFs
784 if (mod_fourier(nloctyp).gt.0) then
787 if (mask_e0(iti,itj).gt.0) then
793 if (mask_ello_PMF.gt.0) then
797 if (mask_turn_PMF.gt.0) then
829 write (iout,*) "Weights:"
831 write (iout,'(a,f10.5)') wname(i),ww(i)
842 if (.not. in_sc_pair_list(iti,j)) then
843 eps(iti,j)=eps_orig(iti,j)
844 & +0.5d0*(x(ii)-eps_orig(iti,iti))
845 eps(j,iti)=eps(iti,j)
859 if (mod_tor .or. mod_fourier(nloctyp).gt.0) then
860 c write (iout,*) "w2x: mod_tor",mod_tor," mod_fourier",
861 c & mod_fourier(nloctyp)
865 do i=-ntortyp+1,ntortyp-1
866 do j=-ntortyp+1,ntortyp-1
867 if (tor_mode.eq.0) then
869 if (mask_tor(0,j,i,iblock).eq.1) then
871 weitor(0,j,i,iblock)=x(ii)
872 else if (mask_tor(0,j,i,iblock).eq.2) then
873 do k=1,nterm(j,i,iblock)
875 v1(k,j,i,iblock)=x(ii)
877 else if (mask_tor(0,j,i,iblock).gt.2) then
878 do k=1,nterm(j,i,iblock)
880 v1(k,j,i,iblock)=x(ii)
882 do k=1,nterm(j,i,iblock)
884 v2(k,j,i,iblock)=x(ii)
889 do k=1,nterm(i,j,iblock)
890 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
891 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
892 v0ij=v0ij+si*v1(k,i,j,iblock)
895 do k=1,nlor(i,j,iblock)
896 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
899 v0(-i,-j,iblock)=v0ij
901 else if (tor_mode.eq.1) then
902 if (mask_tor(0,j,i,1).eq.1) then
904 weitor(0,j,i,1)=x(ii)
905 else if (mask_tor(0,j,i,1).eq.2) then
906 do k=1,nterm_kcc(j,i)
907 do l=1,nterm_kcc_Tb(j,i)
908 do ll=1,nterm_kcc_Tb(j,i)
910 v1_kcc(ll,l,k,j,i)=x(ii)
914 else if (mask_tor(0,j,i,1).gt.2) then
915 do k=1,nterm_kcc(j,i)
916 do l=1,nterm_kcc_Tb(j,i)
917 do ll=1,nterm_kcc_Tb(j,i)
919 v1_kcc(ll,l,k,j,i)=x(ii)
923 do k=1,nterm_kcc(j,i)
924 do l=1,nterm_kcc_Tb(j,i)
925 do ll=1,nterm_kcc_Tb(j,i)
927 v2_kcc(ll,l,k,j,i)=x(ii)
934 c Compute torsional coefficients from local energy surface expansion coeffs
935 if (mask_tor(0,j,i,1).eq.1) then
937 weitor(0,j,i,1)=x(ii)
938 c else if (mask_tor(0,j,i,1).eq.2) then
941 write (iout,*) "ERROR: TOR_MODE>1 ALLOWED ONLY WITH NEWCORR"
951 do i=-nsccortyp,nsccortyp
952 do j=-nsccortyp,nsccortyp
954 if (mask_tor(k,j,i,1).eq.1) then
956 weitor(k,j,i,1)=x(ii)
957 else if (mask_tor(k,j,i,1).eq.2) then
960 v1sccor(l,k,j,i)=x(ii)
962 else if (mask_tor(k,j,i,1).gt.2) then
963 do l=1,nterm_sccor(j,i)
965 v1sccor(l,k,j,i)=x(ii)
967 do l=1,nterm_sccor(j,i)
969 v2sccor(l,k,j,i)=x(ii)
978 c 7/8/17 AL Optimization of the bending potentials
982 if (mask_ang(i).eq.0) cycle
983 do j=1,nbend_kcc_TB(i)
985 v1bend_chyb(j,i)=x(ii)
986 v1bend_chyb(j,-i)=x(ii)
992 if (mod_fourier(nloctyp).gt.0) then
997 if (mod_fourier(i).gt.0) then
1001 if (mask_bnew1(k,j,i).gt.0) then
1010 if (mask_bnew2(k,j,i).gt.0) then
1019 if (mask_ccnew(k,j,i).gt.0) then
1028 if (mask_ddnew(k,j,i).gt.0) then
1036 if (mask_e0new(j,i).gt.0) then
1045 if (mask_eenew(l,k,j,i).gt.0) then
1047 eenew(l,k,j,i)=x(ii)
1059 bnew1(j,1,-i)= bnew1(j,1,i)
1060 bnew1(j,2,-i)=-bnew1(j,2,i)
1061 bnew2(j,1,-i)= bnew2(j,1,i)
1062 bnew2(j,2,-i)=-bnew2(j,2,i)
1065 c ccnew(j,1,i)=ccnew(j,1,i)/2
1066 c ccnew(j,2,i)=ccnew(j,2,i)/2
1067 c ddnew(j,1,i)=ddnew(j,1,i)/2
1068 c ddnew(j,2,i)=ddnew(j,2,i)/2
1069 ccnew(j,1,-i)=ccnew(j,1,i)
1070 ccnew(j,2,-i)=-ccnew(j,2,i)
1071 ddnew(j,1,-i)=ddnew(j,1,i)
1072 ddnew(j,2,-i)=-ddnew(j,2,i)
1074 e0new(1,-i)= e0new(1,i)
1075 e0new(2,-i)=-e0new(2,i)
1076 e0new(3,-i)=-e0new(3,i)
1078 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1079 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1080 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1081 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1085 write (iout,'(a)') "Coefficients of the multibody terms"
1086 do i=-nloctyp+1,nloctyp-1
1087 write (iout,*) "Type: ",onelet(iloctyp(i))
1088 write (iout,*) "Coefficients of the expansion of B1"
1090 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1092 write (iout,*) "Coefficients of the expansion of B2"
1094 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1096 write (iout,*) "Coefficients of the expansion of C"
1097 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1098 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1099 write (iout,*) "Coefficients of the expansion of D"
1100 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1101 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1102 write (iout,*) "Coefficients of the expansion of E"
1103 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1106 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1113 IF (SPLIT_FOURIERTOR .and. tor_mode.eq.2) THEN
1116 if (mod_fouriertor(i).gt.0) then
1120 if (mask_bnew1tor(k,j,i).gt.0) then
1122 bnew1tor(k,j,i)=x(ii)
1129 if (mask_bnew2tor(k,j,i).gt.0) then
1131 bnew2tor(k,j,i)=x(ii)
1138 if (mask_ccnewtor(k,j,i).gt.0) then
1140 ccnewtor(k,j,i)=x(ii)
1147 if (mask_ddnewtor(k,j,i).gt.0) then
1149 ddnewtor(k,j,i)=x(ii)
1155 if (mask_e0newtor(j,i).gt.0) then
1164 if (mask_eenewtor(l,k,j,i).gt.0) then
1166 eenewtor(l,k,j,i)=x(ii)
1179 bnew1tor(j,1,-i)= bnew1tor(j,1,i)
1180 bnew1tor(j,2,-i)=-bnew1tor(j,2,i)
1181 bnew2tor(j,1,-i)= bnew2tor(j,1,i)
1182 bnew2tor(j,2,-i)=-bnew2tor(j,2,i)
1185 c ccnew(j,1,i)=ccnew(j,1,i)/2
1186 c ccnew(j,2,i)=ccnew(j,2,i)/2
1187 c ddnew(j,1,i)=ddnew(j,1,i)/2
1188 c ddnew(j,2,i)=ddnew(j,2,i)/2
1189 ccnewtor(j,1,-i)=ccnewtor(j,1,i)
1190 ccnewtor(j,2,-i)=-ccnewtor(j,2,i)
1191 ddnewtor(j,1,-i)=ddnewtor(j,1,i)
1192 ddnewtor(j,2,-i)=-ddnewtor(j,2,i)
1194 e0newtor(1,-i)= e0newtor(1,i)
1195 e0newtor(2,-i)=-e0newtor(2,i)
1196 e0newtor(3,-i)=-e0newtor(3,i)
1198 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1199 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1200 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1201 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1206 & "Coefficients of the single-body torsional terms"
1207 do i=-nloctyp+1,nloctyp-1
1208 write (iout,*) "Type: ",onelet(iloctyp(i))
1209 write (iout,*) "Coefficients of the expansion of B1"
1211 write (iout,'(3hB1(,i1,1h),3f10.5)')
1212 & j,(bnew1tor(k,j,i),k=1,3)
1214 write (iout,*) "Coefficients of the expansion of B2"
1216 write (iout,'(3hB2(,i1,1h),3f10.5)')
1217 & j,(bnew2tor(k,j,i),k=1,3)
1219 write (iout,*) "Coefficients of the expansion of C"
1220 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1221 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1222 write (iout,*) "Coefficients of the expansion of D"
1223 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1224 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1225 write (iout,*) "Coefficients of the expansion of E"
1226 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1229 write (iout,'(1hE,2i1,2f10.5)')
1230 & j,k,(eenewtor(l,j,k,i),l=1,2)
1239 do i=-nloctyp+1,nloctyp-1
1242 bnew1tor(k,j,i)=bnew1(k,j,i)
1247 bnew2tor(k,j,i)=bnew2(k,j,i)
1251 ccnewtor(k,1,i)=ccnew(k,1,i)
1252 ccnewtor(k,2,i)=ccnew(k,2,i)
1253 ddnewtor(k,1,i)=ddnew(k,1,i)
1254 ddnewtor(k,2,i)=ddnew(k,2,i)
1258 e0newtor(j,i)=e0new(j,i)
1264 eenewtor(l,k,j,i)=eenew(l,k,j,i)
1271 ENDIF ! SPLIT_FOURIERTOR
1273 if (tor_mode.eq.2) then
1274 c Recalculate new torsional potential coefficients
1275 do i=-ntortyp+1,ntortyp-1
1276 do j=-ntortyp+1,ntortyp-1
1277 do k=1,nterm_kcc_Tb(j,i)
1278 do l=1,nterm_kcc_Tb(j,i)
1279 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)
1280 & +bnew1tor(k,2,i)*bnew2tor(l,2,j)
1281 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)
1282 & +bnew1tor(k,2,i)*bnew2tor(l,1,j)
1285 do k=1,nterm_kcc_Tb(j,i)
1286 do l=1,nterm_kcc_Tb(j,i)
1287 c v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1288 c & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1289 c v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1290 c & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1291 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j)
1292 & -ccnewtor(k,2,i)*ddnewtor(l,2,j))
1293 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j)
1294 & +ccnewtor(k,1,i)*ddnewtor(l,2,j))
1303 if (mod_fourier(i).gt.0) then
1306 if (mask_fourier(j,i).gt.0) then
1313 write (iout,*) "Fourier coefficients of cumulants"
1314 write (iout,*) 'Type',i
1315 write (iout,'(a,i2,a,f10.5)') ('b(',j,')=',b(j,i),j=1,13)
1326 c B1tilde(1,i) = b(3,i)
1327 c B1tilde(2,i) =-b(5,i)
1330 CCold(1,1,i)= b(7,i)
1331 CCold(2,2,i)=-b(7,i)
1332 CCold(2,1,i)= b(9,i)
1333 CCold(1,2,i)= b(9,i)
1334 CCold(1,1,-i)= b(7,i)
1335 CCold(2,2,-i)=-b(7,i)
1336 CCold(2,1,-i)=-b(9,i)
1337 CCold(1,2,-i)=-b(9,i)
1338 c Ctilde(1,1,i)=b(7,i)
1339 c Ctilde(1,2,i)=b(9,i)
1340 c Ctilde(2,1,i)=-b(9,i)
1341 c Ctilde(2,2,i)=b(7,i)
1342 DDold(1,1,i)= b(6,i)
1343 DDold(2,2,i)=-b(6,i)
1344 DDold(2,1,i)= b(8,i)
1345 DDold(1,2,i)= b(8,i)
1346 DDold(1,1,-i)= b(6,i)
1347 DDold(2,2,-i)=-b(6,i)
1348 DDold(2,1,-i)=-b(8,i)
1349 DDold(1,2,-i)=-b(8,i)
1350 c Dtilde(1,1,i)=b(6,i)
1351 c Dtilde(1,2,i)=b(8,i)
1352 c Dtilde(2,1,i)=-b(8,i)
1353 c Dtilde(2,2,i)=b(6,i)
1354 EEold(1,1,i)= b(10,i)+b(11,i)
1355 EEold(2,2,i)=-b(10,i)+b(11,i)
1356 EEold(2,1,i)= b(12,i)-b(13,i)
1357 EEold(1,2,i)= b(12,i)+b(13,i)
1358 EEold(1,1,i)= b(10,i)+b(11,i)
1359 EEold(2,2,i)=-b(10,i)+b(11,i)
1360 EEold(2,1,i)= b(12,i)-b(13,i)
1361 EEold(1,2,i)= b(12,i)+b(13,i)
1362 EEold(1,1,-i)= b(10,i)+b(11,i)
1363 EEold(2,2,-i)=-b(10,i)+b(11,i)
1364 EEold(2,1,-i)=-b(12,i)+b(13,i)
1365 EEold(1,2,-i)=-b(12,i)-b(13,i)
1373 &"Coefficients of the cumulants (independent of valence angles)"
1374 do i=-nloctyp+1,nloctyp-1
1375 write (iout,*) 'Type ',onelet(iloctyp(i))
1377 write(iout,'(2f10.5)') B(3,i),B(5,i)
1379 write(iout,'(2f10.5)') B(2,i),B(4,i)
1382 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1386 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1390 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1397 c SC sigmas and anisotropies as parameters
1398 if (mod_side_other) then
1400 c write (iout,*) "mask",i,(mask_sigma(k,i),k=1,4)
1403 if (mask_sigma(1,i).gt.0) then
1406 c write(iout,*) "sig0",i,ii,x(ii)
1408 if (mask_sigma(2,i).gt.0) then
1411 c write(iout,*) "sigi",i,ii,x(ii)
1413 if (mask_sigma(3,i).gt.0) then
1416 c write(iout,*) "chip",i,ii,x(ii)
1418 if (mask_sigma(4,i).gt.0) then
1421 c write(iout,*) "alp",i,ii,x(ii)
1426 if (mod_side.or.mod_side_other) then
1427 if (lprint) write(iout,'(/3a,2i3)') 'Potential is ',
1428 & potname(ipot),', exponents are ',expon,2*expon
1429 C-----------------------------------------------------------------------
1430 C Calculate the "working" parameters of SC interactions.
1431 dwa16=2.0d0**(1.0d0/6.0d0)
1434 chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
1439 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
1440 sigma(j,i)=sigma(i,j)
1441 rs0(i,j)=dwa16*sigma(i,j)
1445 if (lprint) call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1446 if (lprint) write (iout,'(/a/10x,7a/72(1h-))')
1447 & 'Working parameters of the SC interactions:',
1448 & ' a ',' b ',' augm ',' sigma ',' r0 ',
1453 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
1462 c sigeps=dsign(1.0D0,epsij)
1463 aa(i,j)=ftune_eps(epsij)*rrij*rrij
1468 sigt1sq=sigma0(i)**2
1469 sigt2sq=sigma0(j)**2
1472 ratsig1=sigt2sq/sigt1sq
1473 ratsig2=1.0D0/ratsig1
1474 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
1475 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
1476 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
1480 sigmaii(i,j)=rsum_max
1481 sigmaii(j,i)=rsum_max
1482 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max)
1484 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
1485 augm(i,j)=epsij*r_augm**(2*expon)
1486 c augm(i,j)=0.5D0**(2*expon)*aa(i,j)
1493 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
1494 & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
1495 & sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
1506 if (mask_elec(i,j,1).gt.0) then
1509 if (epp(i,j).lt.0.0d0) epp(i,j)=0.0d0
1512 if (mask_elec(i,j,2).gt.0) then
1514 rpp(i,j)=dabs(x(ii))
1515 rpp(j,i)=dabs(x(ii))
1517 if (mask_elec(i,j,3).gt.0) then
1520 if (elpp6(i,j).gt.0.0d0) elpp6(i,j)=0.0d0
1521 elpp6(j,i)=elpp6(i,j)
1523 if (mask_elec(i,j,4).gt.0) then
1529 app (i,j)=epp(i,j)*rri*rri
1530 app (j,i)=epp(j,i)*rri*rri
1531 bpp (i,j)=-2.0D0*epp(i,j)*rri
1532 bpp (j,i)=-2.0D0*epp(i,j)*rri
1533 ael6(i,j)=elpp6(i,j)*4.2D0**6
1534 ael6(j,i)=elpp6(i,j)*4.2D0**6
1535 ael3(i,j)=elpp3(i,j)*4.2D0**3
1536 ael3(j,i)=elpp3(i,j)*4.2D0**3
1540 write (iout,'(/a)') 'Electrostatic interaction constants:'
1541 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)')
1542 & 'IT','JT','APP','BPP','AEL6','AEL3'
1545 write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
1546 & ael6(i,j),ael3(i,j)
1553 C Define the SC-p interaction constants
1557 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1558 & mask_scp(0,2,2) .gt. 0) then
1560 if (mask_scp(0,j,1).gt.0) then
1566 if (mask_scp(0,j,2).gt.0) then
1576 if (mask_scp(i,j,1).gt.0) then
1580 if (mask_scp(i,j,2).gt.0) then
1588 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
1589 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
1590 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
1591 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
1595 write (iout,*) "Parameters of SC-p interactions:"
1597 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1598 & eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1604 c Weight of the shielding term
1605 if (imask(ind_shield).eq.1) then
1607 ww(ind_shield)=x(ii)
1609 wshield=ww(ind_shield)
1611 c write (iout,*) "x2w: nvar",nvarr
1613 c Base of heat capacity
1616 if (nbeta(2,iprot).gt.0) then
1618 heatbase(iprot)=x(ii)
1622 if (ii.ne.nvarr) then
1623 write (iout,*) "!!!!! SEVERE ERROR IN X2W !!!!! (chuj nastapil)"
1624 write (iout,*) "NVAR=",nvarr," LAST VARIABLE INDEX=",ii
1625 stop "CHUJ ABSOLUTNY!!!!"
1630 c---------------------------------------------------------------------------
1631 double precision function ftune_eps(x)
1633 double precision x,y
1634 double precision delta /1.0d-5/
1635 if (dabs(x).lt.delta) then
1637 ftune_eps=delta*(3.0d0+6.0d0*y**2-y**4)/8
1643 c---------------------------------------------------------------------------
1644 double precision function ftune_epsprim(x)
1646 double precision x,y
1647 double precision delta /1.0d-5/
1648 if (dabs(x).lt.delta) then
1650 ftune_epsprim=1.5d0*y-0.5*y**3
1652 ftune_epsprim=dsign(1.0d0,x)
1656 c---------------------------------------------------------------------------
1657 double precision function ftune_epsbis(x)
1659 double precision x,y
1660 double precision delta /1.0d-5/
1661 if (dabs(x).lt.delta) then
1663 ftune_epsbis=(1.5d0-1.5d0*y**3)/delta
1669 c---------------------------------------------------------------------------
1670 logical function in_sc_pair_list(iti,itj)
1672 include "DIMENSIONS"
1673 include "DIMENSIONS.ZSCOPT"
1674 include "COMMON.WEIGHTS"
1677 if (ityp_psc(1,i).eq.iti .and. ityp_psc(2,i).eq.itj .or.
1678 & ityp_psc(1,i).eq.itj .and. ityp_psc(2,i).eq.iti) then
1679 in_sc_pair_list=.true.
1683 in_sc_pair_list=.false.