8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
873 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
874 if (oldion.eq.1) then
876 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
877 print *,msc(i,5),restok(i,5)
881 read (iion,*) ncatprotparm
882 allocate(catprm(ncatprotparm,4))
884 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
888 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
892 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
895 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
896 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
900 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
901 restyp(i,5),"-",restyp(j,2)
904 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
905 !----------------------------------------------------
906 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
907 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
908 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
909 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
910 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
911 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
921 allocate(liptranene(ntyp))
922 !C reading lipid parameters
923 write (iout,*) "iliptranpar",iliptranpar
925 read(iliptranpar,*) pepliptran
928 read(iliptranpar,*) liptranene(i)
929 print *,liptranene(i)
935 ! Read the parameters of the probability distribution/energy expression
936 ! of the virtual-bond valence angles theta
939 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
940 (bthet(j,i,1,1),j=1,2)
941 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
942 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
943 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
947 athet(1,i,1,-1)=athet(1,i,1,1)
948 athet(2,i,1,-1)=athet(2,i,1,1)
949 bthet(1,i,1,-1)=-bthet(1,i,1,1)
950 bthet(2,i,1,-1)=-bthet(2,i,1,1)
951 athet(1,i,-1,1)=-athet(1,i,1,1)
952 athet(2,i,-1,1)=-athet(2,i,1,1)
953 bthet(1,i,-1,1)=bthet(1,i,1,1)
954 bthet(2,i,-1,1)=bthet(2,i,1,1)
958 athet(1,i,-1,-1)=athet(1,-i,1,1)
959 athet(2,i,-1,-1)=-athet(2,-i,1,1)
960 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
961 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
962 athet(1,i,-1,1)=athet(1,-i,1,1)
963 athet(2,i,-1,1)=-athet(2,-i,1,1)
964 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
965 bthet(2,i,-1,1)=bthet(2,-i,1,1)
966 athet(1,i,1,-1)=-athet(1,-i,1,1)
967 athet(2,i,1,-1)=athet(2,-i,1,1)
968 bthet(1,i,1,-1)=bthet(1,-i,1,1)
969 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
974 polthet(j,i)=polthet(j,-i)
977 gthet(j,i)=gthet(j,-i)
985 'Parameters of the virtual-bond valence angles:'
986 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
987 ' ATHETA0 ',' A1 ',' A2 ',&
990 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
991 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
993 write (iout,'(/a/9x,5a/79(1h-))') &
994 'Parameters of the expression for sigma(theta_c):',&
995 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
996 ' ALPH3 ',' SIGMA0C '
998 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
999 (polthet(j,i),j=0,3),sigc0(i)
1001 write (iout,'(/a/9x,5a/79(1h-))') &
1002 'Parameters of the second gaussian:',&
1003 ' THETA0 ',' SIGMA0 ',' G1 ',&
1006 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1007 sig0(i),(gthet(j,i),j=1,3)
1010 write (iout,'(a)') &
1011 'Parameters of the virtual-bond valence angles:'
1012 write (iout,'(/a/9x,5a/79(1h-))') &
1013 'Coefficients of expansion',&
1014 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1015 ' b1*10^1 ',' b2*10^1 '
1017 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1018 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1019 (10*bthet(j,i,1,1),j=1,2)
1021 write (iout,'(/a/9x,5a/79(1h-))') &
1022 'Parameters of the expression for sigma(theta_c):',&
1023 ' alpha0 ',' alph1 ',' alph2 ',&
1024 ' alhp3 ',' sigma0c '
1026 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1027 (polthet(j,i),j=0,3),sigc0(i)
1029 write (iout,'(/a/9x,5a/79(1h-))') &
1030 'Parameters of the second gaussian:',&
1031 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1034 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1035 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1041 ! Read the parameters of Utheta determined from ab initio surfaces
1042 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1044 IF (tor_mode.eq.0) THEN
1045 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1046 ntheterm3,nsingle,ndouble
1047 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1049 !----------------------------------------------------
1050 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1051 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1052 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1053 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1054 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1055 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1056 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1057 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1058 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1059 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1060 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1061 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1062 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1063 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1064 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1065 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1066 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1067 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1068 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1069 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1070 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1071 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1072 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1073 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1075 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1077 ithetyp(i)=-ithetyp(-i)
1080 aa0thet(:,:,:,:)=0.0d0
1081 aathet(:,:,:,:,:)=0.0d0
1082 bbthet(:,:,:,:,:,:)=0.0d0
1083 ccthet(:,:,:,:,:,:)=0.0d0
1084 ddthet(:,:,:,:,:,:)=0.0d0
1085 eethet(:,:,:,:,:,:)=0.0d0
1086 ffthet(:,:,:,:,:,:,:)=0.0d0
1087 ggthet(:,:,:,:,:,:,:)=0.0d0
1089 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1091 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1092 ! VAR:1=non-glicyne non-proline 2=proline
1093 ! VAR:negative values for D-aminoacid
1095 do j=-nthetyp,nthetyp
1096 do k=-nthetyp,nthetyp
1097 read (ithep,'(6a)',end=111,err=111) res1
1098 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1099 ! VAR: aa0thet is variable describing the average value of Foureir
1100 ! VAR: expansion series
1101 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1102 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1103 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1104 read (ithep,*,end=111,err=111) &
1105 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1106 read (ithep,*,end=111,err=111) &
1107 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1108 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1109 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1110 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1112 read (ithep,*,end=111,err=111) &
1113 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1114 ffthet(lll,llll,ll,i,j,k,iblock),&
1115 ggthet(llll,lll,ll,i,j,k,iblock),&
1116 ggthet(lll,llll,ll,i,j,k,iblock),&
1117 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1122 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1123 ! coefficients of theta-and-gamma-dependent terms are zero.
1124 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1125 ! RECOMENTDED AFTER VERSION 3.3)
1129 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1130 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1132 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1133 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1136 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1138 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1141 ! AND COMMENT THE LOOPS BELOW
1145 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1146 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1148 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1149 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1152 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1154 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1159 ! Substitution for D aminoacids from symmetry.
1162 do j=-nthetyp,nthetyp
1163 do k=-nthetyp,nthetyp
1164 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1166 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1170 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1171 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1172 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1173 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1179 ffthet(llll,lll,ll,i,j,k,iblock)= &
1180 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1181 ffthet(lll,llll,ll,i,j,k,iblock)= &
1182 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1183 ggthet(llll,lll,ll,i,j,k,iblock)= &
1184 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1185 ggthet(lll,llll,ll,i,j,k,iblock)= &
1186 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1195 ! Control printout of the coefficients of virtual-bond-angle potentials
1198 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1203 write (iout,'(//4a)') &
1204 'Type ',onelett(i),onelett(j),onelett(k)
1205 write (iout,'(//a,10x,a)') " l","a[l]"
1206 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1207 write (iout,'(i2,1pe15.5)') &
1208 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1210 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1211 "b",l,"c",l,"d",l,"e",l
1213 write (iout,'(i2,4(1pe15.5))') m,&
1214 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1215 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1219 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1220 "f+",l,"f-",l,"g+",l,"g-",l
1223 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1224 ffthet(n,m,l,i,j,k,iblock),&
1225 ffthet(m,n,l,i,j,k,iblock),&
1226 ggthet(n,m,l,i,j,k,iblock),&
1227 ggthet(m,n,l,i,j,k,iblock)
1238 !C here will be the apropriate recalibrating for D-aminoacid
1239 read (ithep,*,end=121,err=121) nthetyp
1240 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1241 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1242 do i=-nthetyp+1,nthetyp-1
1243 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1244 do j=0,nbend_kcc_Tb(i)
1245 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1249 write (iout,'(a)') &
1250 "Parameters of the valence-only potentials"
1251 do i=-nthetyp+1,nthetyp-1
1252 write (iout,'(2a)') "Type ",toronelet(i)
1253 do k=0,nbend_kcc_Tb(i)
1254 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1260 write (2,*) "Start reading THETA_PDB",ithep_pdb
1262 ! write (2,*) 'i=',i
1263 read (ithep_pdb,*,err=111,end=111) &
1264 a0thet(i),(athet(j,i,1,1),j=1,2),&
1265 (bthet(j,i,1,1),j=1,2)
1266 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1267 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1268 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1269 sigc0(i)=sigc0(i)**2
1272 athet(1,i,1,-1)=athet(1,i,1,1)
1273 athet(2,i,1,-1)=athet(2,i,1,1)
1274 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1275 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1276 athet(1,i,-1,1)=-athet(1,i,1,1)
1277 athet(2,i,-1,1)=-athet(2,i,1,1)
1278 bthet(1,i,-1,1)=bthet(1,i,1,1)
1279 bthet(2,i,-1,1)=bthet(2,i,1,1)
1282 a0thet(i)=a0thet(-i)
1283 athet(1,i,-1,-1)=athet(1,-i,1,1)
1284 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1285 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1286 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1287 athet(1,i,-1,1)=athet(1,-i,1,1)
1288 athet(2,i,-1,1)=-athet(2,-i,1,1)
1289 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1290 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1291 athet(1,i,1,-1)=-athet(1,-i,1,1)
1292 athet(2,i,1,-1)=athet(2,-i,1,1)
1293 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1294 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1295 theta0(i)=theta0(-i)
1299 polthet(j,i)=polthet(j,-i)
1302 gthet(j,i)=gthet(j,-i)
1305 write (2,*) "End reading THETA_PDB"
1309 !--------------- Reading theta parameters for nucleic acid-------
1310 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1311 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1312 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1313 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1314 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1315 nthetyp_nucl+1,nthetyp_nucl+1))
1316 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1317 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1318 nthetyp_nucl+1,nthetyp_nucl+1))
1319 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1320 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1321 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1322 nthetyp_nucl+1,nthetyp_nucl+1))
1323 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1324 nthetyp_nucl+1,nthetyp_nucl+1))
1325 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1326 nthetyp_nucl+1,nthetyp_nucl+1))
1327 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1328 nthetyp_nucl+1,nthetyp_nucl+1))
1329 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1330 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1331 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1332 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1333 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1334 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1336 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1337 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1339 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1341 aa0thet_nucl(:,:,:)=0.0d0
1342 aathet_nucl(:,:,:,:)=0.0d0
1343 bbthet_nucl(:,:,:,:,:)=0.0d0
1344 ccthet_nucl(:,:,:,:,:)=0.0d0
1345 ddthet_nucl(:,:,:,:,:)=0.0d0
1346 eethet_nucl(:,:,:,:,:)=0.0d0
1347 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1348 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1353 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1354 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1355 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1356 read (ithep_nucl,*,end=111,err=111) &
1357 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1358 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1359 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1360 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1361 read (ithep_nucl,*,end=111,err=111) &
1362 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1363 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1364 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1369 !-------------------------------------------
1370 allocate(nlob(ntyp1)) !(ntyp1)
1371 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1372 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1373 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1380 gaussc(:,:,:,:)=0.0D0
1384 ! Read the parameters of the probability distribution/energy expression
1385 ! of the side chains.
1388 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1392 dsc_inv(i)=1.0D0/dsc(i)
1403 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1404 ((blower(k,l,1),l=1,k),k=1,3)
1405 censc(1,1,-i)=censc(1,1,i)
1406 censc(2,1,-i)=censc(2,1,i)
1407 censc(3,1,-i)=-censc(3,1,i)
1409 read (irotam,*,end=112,err=112) bsc(j,i)
1410 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1411 ((blower(k,l,j),l=1,k),k=1,3)
1412 censc(1,j,-i)=censc(1,j,i)
1413 censc(2,j,-i)=censc(2,j,i)
1414 censc(3,j,-i)=-censc(3,j,i)
1415 ! BSC is amplitude of Gaussian
1422 akl=akl+blower(k,m,j)*blower(l,m,j)
1426 if (((k.eq.3).and.(l.ne.3)) &
1427 .or.((l.eq.3).and.(k.ne.3))) then
1428 gaussc(k,l,j,-i)=-akl
1429 gaussc(l,k,j,-i)=-akl
1431 gaussc(k,l,j,-i)=akl
1432 gaussc(l,k,j,-i)=akl
1441 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1444 if (nlobi.gt.0) then
1446 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1447 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1448 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1449 'log h',(bsc(j,i),j=1,nlobi)
1450 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1451 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1453 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1454 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1457 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1458 write (iout,'(a,f10.4,4(16x,f10.4))') &
1459 'Center ',(bsc(j,i),j=1,nlobi)
1460 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1469 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1470 ! added by Urszula Kozlowska 07/11/2007
1472 !el Maximum number of SC local term fitting function coefficiants
1473 !el integer,parameter :: maxsccoef=65
1475 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1478 read (irotam,*,end=112,err=112)
1480 read (irotam,*,end=112,err=112)
1483 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1487 !---------reading nucleic acid parameters for rotamers-------------------
1488 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1489 do i=1,ntyp_molec(2)
1490 read (irotam_nucl,*,end=112,err=112)
1492 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1498 write (iout,*) "Base rotamer parameters"
1499 do i=1,ntyp_molec(2)
1500 write (iout,'(a)') restyp(i,2)
1501 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1506 ! Read the parameters of the probability distribution/energy expression
1507 ! of the side chains.
1509 write (2,*) "Start reading ROTAM_PDB"
1511 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1515 dsc_inv(i)=1.0D0/dsc(i)
1526 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1527 ((blower(k,l,1),l=1,k),k=1,3)
1529 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1530 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1531 ((blower(k,l,j),l=1,k),k=1,3)
1538 akl=akl+blower(k,m,j)*blower(l,m,j)
1548 write (2,*) "End reading ROTAM_PDB"
1554 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1555 !C interaction energy of the Gly, Ala, and Pro prototypes.
1557 read (ifourier,*) nloctyp
1558 SPLIT_FOURIERTOR = nloctyp.lt.0
1559 nloctyp = iabs(nloctyp)
1560 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1561 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1562 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1563 !C allocate(ctilde(2,2,nres))
1564 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1565 !C allocate(gtb1(2,nres))
1566 !C allocate(gtb2(2,nres))
1567 !C allocate(cc(2,2,nres))
1568 !C allocate(dd(2,2,nres))
1569 !C allocate(ee(2,2,nres))
1570 !C allocate(gtcc(2,2,nres))
1571 !C allocate(gtdd(2,2,nres))
1572 !C allocate(gtee(2,2,nres))
1575 allocate(itype2loc(-ntyp1:ntyp1))
1576 allocate(iloctyp(-nloctyp:nloctyp))
1577 allocate(bnew1(3,2,-nloctyp:nloctyp))
1578 allocate(bnew2(3,2,-nloctyp:nloctyp))
1579 allocate(ccnew(3,2,-nloctyp:nloctyp))
1580 allocate(ddnew(3,2,-nloctyp:nloctyp))
1581 allocate(e0new(3,-nloctyp:nloctyp))
1582 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1583 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1584 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1585 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1586 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1587 allocate(e0newtor(3,-nloctyp:nloctyp))
1588 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1590 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1591 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1592 itype2loc(ntyp1)=nloctyp
1593 iloctyp(nloctyp)=ntyp1
1595 itype2loc(-i)=-itype2loc(i)
1598 allocate(iloctyp(-nloctyp:nloctyp))
1599 allocate(itype2loc(-ntyp1:ntyp1))
1606 iloctyp(-i)=-iloctyp(i)
1608 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1609 !c write (iout,*) "nloctyp",nloctyp,
1610 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1611 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1612 !c write (iout,*) "nloctyp",nloctyp,
1613 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1616 !c write (iout,*) "NEWCORR",i
1617 read (ifourier,*,end=115,err=115)
1620 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1623 !c write (iout,*) "NEWCORR BNEW1"
1624 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1627 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1630 !c write (iout,*) "NEWCORR BNEW2"
1631 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1633 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1634 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1636 !c write (iout,*) "NEWCORR CCNEW"
1637 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1639 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1640 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1642 !c write (iout,*) "NEWCORR DDNEW"
1643 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1647 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1651 !c write (iout,*) "NEWCORR EENEW1"
1652 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1654 read (ifourier,*,end=115,err=115) e0new(ii,i)
1656 !c write (iout,*) (e0new(ii,i),ii=1,3)
1658 !c write (iout,*) "NEWCORR EENEW"
1661 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1662 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1663 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1664 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1669 bnew1(ii,1,-i)= bnew1(ii,1,i)
1670 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1671 bnew2(ii,1,-i)= bnew2(ii,1,i)
1672 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1675 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1676 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1677 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1678 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1679 ccnew(ii,1,-i)=ccnew(ii,1,i)
1680 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1681 ddnew(ii,1,-i)=ddnew(ii,1,i)
1682 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1684 e0new(1,-i)= e0new(1,i)
1685 e0new(2,-i)=-e0new(2,i)
1686 e0new(3,-i)=-e0new(3,i)
1688 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1689 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1690 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1691 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1695 write (iout,'(a)') "Coefficients of the multibody terms"
1696 do i=-nloctyp+1,nloctyp-1
1697 write (iout,*) "Type: ",onelet(iloctyp(i))
1698 write (iout,*) "Coefficients of the expansion of B1"
1700 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1702 write (iout,*) "Coefficients of the expansion of B2"
1704 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1706 write (iout,*) "Coefficients of the expansion of C"
1707 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1708 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1709 write (iout,*) "Coefficients of the expansion of D"
1710 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1711 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1712 write (iout,*) "Coefficients of the expansion of E"
1713 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1716 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1721 IF (SPLIT_FOURIERTOR) THEN
1723 !c write (iout,*) "NEWCORR TOR",i
1724 read (ifourier,*,end=115,err=115)
1727 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1730 !c write (iout,*) "NEWCORR BNEW1 TOR"
1731 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1734 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1737 !c write (iout,*) "NEWCORR BNEW2 TOR"
1738 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1740 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1741 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1743 !c write (iout,*) "NEWCORR CCNEW TOR"
1744 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1746 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1747 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1749 !c write (iout,*) "NEWCORR DDNEW TOR"
1750 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1754 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1758 !c write (iout,*) "NEWCORR EENEW1 TOR"
1759 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1761 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1763 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1765 !c write (iout,*) "NEWCORR EENEW TOR"
1768 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1769 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1770 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1771 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1776 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1777 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1778 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1779 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1782 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1783 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1784 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1785 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1787 e0newtor(1,-i)= e0newtor(1,i)
1788 e0newtor(2,-i)=-e0newtor(2,i)
1789 e0newtor(3,-i)=-e0newtor(3,i)
1791 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1792 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1793 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1794 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1798 write (iout,'(a)') &
1799 "Single-body coefficients of the torsional potentials"
1800 do i=-nloctyp+1,nloctyp-1
1801 write (iout,*) "Type: ",onelet(iloctyp(i))
1802 write (iout,*) "Coefficients of the expansion of B1tor"
1804 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1805 j,(bnew1tor(k,j,i),k=1,3)
1807 write (iout,*) "Coefficients of the expansion of B2tor"
1809 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1810 j,(bnew2tor(k,j,i),k=1,3)
1812 write (iout,*) "Coefficients of the expansion of Ctor"
1813 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1814 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1815 write (iout,*) "Coefficients of the expansion of Dtor"
1816 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1817 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1818 write (iout,*) "Coefficients of the expansion of Etor"
1819 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1822 write (iout,'(1hE,2i1,2f10.5)') &
1823 j,k,(eenewtor(l,j,k,i),l=1,2)
1829 do i=-nloctyp+1,nloctyp-1
1832 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1837 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1841 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1842 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1843 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1844 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1847 ENDIF !SPLIT_FOURIER_TOR
1849 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1850 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1851 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1852 allocate(b(13,-nloctyp-1:nloctyp+1))
1854 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1856 read (ifourier,*,end=115,err=115)
1857 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1859 write (iout,*) 'Type ',onelet(iloctyp(i))
1860 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1874 !c B1tilde(1,i) = b(3)
1875 !c! B1tilde(2,i) =-b(5)
1876 !c! B1tilde(1,-i) =-b(3)
1877 !c! B1tilde(2,-i) =b(5)
1878 !c! b1tilde(1,i)=0.0d0
1879 !c b1tilde(2,i)=0.0d0
1884 !cc B1tilde(1,i) = b(3,i)
1885 !cc B1tilde(2,i) =-b(5,i)
1886 !c B1tilde(1,-i) =-b(3,i)
1887 !c B1tilde(2,-i) =b(5,i)
1888 !cc b1tilde(1,i)=0.0d0
1889 !cc b1tilde(2,i)=0.0d0
1890 !cc B2(1,i) = b(2,i)
1891 !cc B2(2,i) = b(4,i)
1893 !c B2(2,-i) =-b(4,i)
1897 CCold(1,1,i)= b(7,i)
1898 CCold(2,2,i)=-b(7,i)
1899 CCold(2,1,i)= b(9,i)
1900 CCold(1,2,i)= b(9,i)
1901 CCold(1,1,-i)= b(7,i)
1902 CCold(2,2,-i)=-b(7,i)
1903 CCold(2,1,-i)=-b(9,i)
1904 CCold(1,2,-i)=-b(9,i)
1909 !c Ctilde(1,1,i)= CCold(1,1,i)
1910 !c Ctilde(1,2,i)= CCold(1,2,i)
1911 !c Ctilde(2,1,i)=-CCold(2,1,i)
1912 !c Ctilde(2,2,i)=-CCold(2,2,i)
1917 !c Ctilde(1,1,i)= CCold(1,1,i)
1918 !c Ctilde(1,2,i)= CCold(1,2,i)
1919 !c Ctilde(2,1,i)=-CCold(2,1,i)
1920 !c Ctilde(2,2,i)=-CCold(2,2,i)
1922 !c Ctilde(1,1,i)=0.0d0
1923 !c Ctilde(1,2,i)=0.0d0
1924 !c Ctilde(2,1,i)=0.0d0
1925 !c Ctilde(2,2,i)=0.0d0
1926 DDold(1,1,i)= b(6,i)
1927 DDold(2,2,i)=-b(6,i)
1928 DDold(2,1,i)= b(8,i)
1929 DDold(1,2,i)= b(8,i)
1930 DDold(1,1,-i)= b(6,i)
1931 DDold(2,2,-i)=-b(6,i)
1932 DDold(2,1,-i)=-b(8,i)
1933 DDold(1,2,-i)=-b(8,i)
1938 !c Dtilde(1,1,i)= DD(1,1,i)
1939 !c Dtilde(1,2,i)= DD(1,2,i)
1940 !c Dtilde(2,1,i)=-DD(2,1,i)
1941 !c Dtilde(2,2,i)=-DD(2,2,i)
1943 !c Dtilde(1,1,i)=0.0d0
1944 !c Dtilde(1,2,i)=0.0d0
1945 !c Dtilde(2,1,i)=0.0d0
1946 !c Dtilde(2,2,i)=0.0d0
1947 EEold(1,1,i)= b(10,i)+b(11,i)
1948 EEold(2,2,i)=-b(10,i)+b(11,i)
1949 EEold(2,1,i)= b(12,i)-b(13,i)
1950 EEold(1,2,i)= b(12,i)+b(13,i)
1951 EEold(1,1,-i)= b(10,i)+b(11,i)
1952 EEold(2,2,-i)=-b(10,i)+b(11,i)
1953 EEold(2,1,-i)=-b(12,i)+b(13,i)
1954 EEold(1,2,-i)=-b(12,i)-b(13,i)
1955 write(iout,*) "TU DOCHODZE"
1961 !c ee(2,1,i)=ee(1,2,i)
1966 "Coefficients of the cumulants (independent of valence angles)"
1967 do i=-nloctyp+1,nloctyp-1
1968 write (iout,*) 'Type ',onelet(iloctyp(i))
1970 write(iout,'(2f10.5)') B(3,i),B(5,i)
1972 write(iout,'(2f10.5)') B(2,i),B(4,i)
1975 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1979 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1983 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1992 ! Read torsional parameters in old format
1994 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1996 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1997 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1998 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2000 !el from energy module--------
2001 allocate(v1(nterm_old,ntortyp,ntortyp))
2002 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2003 !el---------------------------
2008 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2014 write (iout,'(/a/)') 'Torsional constants:'
2017 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2018 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2024 ! Read torsional parameters
2026 IF (TOR_MODE.eq.0) THEN
2027 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2028 read (itorp,*,end=113,err=113) ntortyp
2029 !el from energy module---------
2030 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2031 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2033 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2034 allocate(vlor2(maxlor,ntortyp,ntortyp))
2035 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2036 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2038 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2039 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2040 !el---------------------------
2043 !el---------------------------
2045 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2047 itortyp(i)=-itortyp(-i)
2049 itortyp(ntyp1)=ntortyp
2050 itortyp(-ntyp1)=-ntortyp
2052 write (iout,*) 'ntortyp',ntortyp
2054 do j=-ntortyp+1,ntortyp-1
2055 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2057 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2058 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2061 do k=1,nterm(i,j,iblock)
2062 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2064 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2065 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2066 v0ij=v0ij+si*v1(k,i,j,iblock)
2068 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2069 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2070 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2072 do k=1,nlor(i,j,iblock)
2073 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2074 vlor2(k,i,j),vlor3(k,i,j)
2075 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2078 v0(-i,-j,iblock)=v0ij
2084 write (iout,'(/a/)') 'Torsional constants:'
2086 do i=-ntortyp,ntortyp
2087 do j=-ntortyp,ntortyp
2088 write (iout,*) 'ityp',i,' jtyp',j
2089 write (iout,*) 'Fourier constants'
2090 do k=1,nterm(i,j,iblock)
2091 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2094 write (iout,*) 'Lorenz constants'
2095 do k=1,nlor(i,j,iblock)
2096 write (iout,'(3(1pe15.5))') &
2097 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2103 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2105 ! 6/23/01 Read parameters for double torsionals
2107 !el from energy module------------
2108 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2109 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2110 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2112 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2113 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2114 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2115 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2116 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2117 !---------------------------------
2121 do j=-ntortyp+1,ntortyp-1
2122 do k=-ntortyp+1,ntortyp-1
2123 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2124 ! write (iout,*) "OK onelett",
2127 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2128 .or. t3.ne.toronelet(k)) then
2129 write (iout,*) "Error in double torsional parameter file",&
2132 call MPI_Finalize(Ierror)
2134 stop "Error in double torsional parameter file"
2136 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2137 ntermd_2(i,j,k,iblock)
2138 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2139 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2140 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2141 ntermd_1(i,j,k,iblock))
2142 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2143 ntermd_1(i,j,k,iblock))
2144 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2145 ntermd_1(i,j,k,iblock))
2146 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2147 ntermd_1(i,j,k,iblock))
2148 ! Martix of D parameters for one dimesional foureir series
2149 do l=1,ntermd_1(i,j,k,iblock)
2150 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2151 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2152 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2153 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2154 ! write(iout,*) "whcodze" ,
2155 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2157 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2158 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2159 v2s(m,l,i,j,k,iblock),&
2160 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2161 ! Martix of D parameters for two dimesional fourier series
2162 do l=1,ntermd_2(i,j,k,iblock)
2164 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2165 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2166 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2167 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2176 write (iout,*) 'Constants for double torsionals'
2179 do j=-ntortyp+1,ntortyp-1
2180 do k=-ntortyp+1,ntortyp-1
2181 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2182 ' nsingle',ntermd_1(i,j,k,iblock),&
2183 ' ndouble',ntermd_2(i,j,k,iblock)
2185 write (iout,*) 'Single angles:'
2186 do l=1,ntermd_1(i,j,k,iblock)
2187 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2188 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2189 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2190 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2193 write (iout,*) 'Pairs of angles:'
2194 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2195 do l=1,ntermd_2(i,j,k,iblock)
2196 write (iout,'(i5,20f10.5)') &
2197 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2200 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2201 do l=1,ntermd_2(i,j,k,iblock)
2202 write (iout,'(i5,20f10.5)') &
2203 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2204 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2214 itype2loc(i)=itortyp(i)
2218 ELSE IF (TOR_MODE.eq.1) THEN
2220 !C read valence-torsional parameters
2221 read (itorp,*,end=121,err=121) ntortyp
2223 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2225 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2226 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2229 itype2loc(i)=itortyp(i)
2233 itortyp(i)=-itortyp(-i)
2235 do i=-ntortyp+1,ntortyp-1
2236 do j=-ntortyp+1,ntortyp-1
2237 !C first we read the cos and sin gamma parameters
2238 read (itorp,'(13x,a)',end=121,err=121) string
2239 write (iout,*) i,j,string
2240 read (itorp,*,end=121,err=121) &
2241 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2242 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2243 do k=1,nterm_kcc(j,i)
2244 do l=1,nterm_kcc_Tb(j,i)
2245 do ll=1,nterm_kcc_Tb(j,i)
2246 read (itorp,*,end=121,err=121) ii,jj,kk, &
2247 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2255 !c AL 4/8/16: Calculate coefficients from one-body parameters
2257 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2258 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2259 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2260 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2261 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2264 print *,i,itortyp(i)
2265 itortyp(i)=itype2loc(i)
2268 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2270 do i=-ntortyp+1,ntortyp-1
2271 do j=-ntortyp+1,ntortyp-1
2274 do k=1,nterm_kcc_Tb(j,i)
2275 do l=1,nterm_kcc_Tb(j,i)
2276 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2277 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2278 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2279 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2282 do k=1,nterm_kcc_Tb(j,i)
2283 do l=1,nterm_kcc_Tb(j,i)
2285 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2286 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2287 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2288 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2290 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2291 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2292 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2293 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2297 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2301 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2306 if (tor_mode.gt.0 .and. lprint) then
2307 !c Print valence-torsional parameters
2308 write (iout,'(a)') &
2309 "Parameters of the valence-torsional potentials"
2310 do i=-ntortyp+1,ntortyp-1
2311 do j=-ntortyp+1,ntortyp-1
2312 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2313 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2314 do k=1,nterm_kcc(j,i)
2315 do l=1,nterm_kcc_Tb(j,i)
2316 do ll=1,nterm_kcc_Tb(j,i)
2317 write (iout,'(3i5,2f15.4)')&
2318 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2326 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2327 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2328 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2329 !el from energy module---------
2330 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2331 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2333 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2334 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2335 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2336 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2338 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2339 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2340 !el---------------------------
2343 !el--------------------
2344 read (itorp_nucl,*,end=113,err=113) &
2345 (itortyp_nucl(i),i=1,ntyp_molec(2))
2346 ! print *,itortyp_nucl(:)
2347 !c write (iout,*) 'ntortyp',ntortyp
2350 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2351 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2354 do k=1,nterm_nucl(i,j)
2355 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2356 v0ij=v0ij+si*v1_nucl(k,i,j)
2359 do k=1,nlor_nucl(i,j)
2360 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2361 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2362 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2368 ! Read of Side-chain backbone correlation parameters
2369 ! Modified 11 May 2012 by Adasko
2372 read (isccor,*,end=119,err=119) nsccortyp
2374 !el from module energy-------------
2375 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2376 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2377 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2378 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2379 !-----------------------------------
2381 !el from module energy-------------
2382 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2384 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2386 isccortyp(i)=-isccortyp(-i)
2388 iscprol=isccortyp(20)
2389 ! write (iout,*) 'ntortyp',ntortyp
2391 !c maxinter is maximum interaction sites
2392 !el from module energy---------
2393 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2394 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2395 -nsccortyp:nsccortyp))
2396 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2397 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2398 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2399 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2400 !-----------------------------------
2404 read (isccor,*,end=119,err=119) &
2405 nterm_sccor(i,j),nlor_sccor(i,j)
2411 nterm_sccor(-i,j)=nterm_sccor(i,j)
2412 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2413 nterm_sccor(i,-j)=nterm_sccor(i,j)
2414 do k=1,nterm_sccor(i,j)
2415 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2417 if (j.eq.iscprol) then
2418 if (i.eq.isccortyp(10)) then
2419 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2420 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2422 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2423 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2424 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2425 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2426 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2427 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2428 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2429 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2432 if (i.eq.isccortyp(10)) then
2433 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2434 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2436 if (j.eq.isccortyp(10)) then
2437 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2438 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2440 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2441 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2442 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2443 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2444 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2445 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2449 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2450 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2451 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2452 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2455 do k=1,nlor_sccor(i,j)
2456 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2457 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2458 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2459 (1+vlor3sccor(k,i,j)**2)
2461 v0sccor(l,i,j)=v0ijsccor
2462 v0sccor(l,-i,j)=v0ijsccor1
2463 v0sccor(l,i,-j)=v0ijsccor2
2464 v0sccor(l,-i,-j)=v0ijsccor3
2470 !el from module energy-------------
2471 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2473 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2474 ! write (iout,*) 'ntortyp',ntortyp
2476 !c maxinter is maximum interaction sites
2477 !el from module energy---------
2478 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2479 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2480 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2481 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2482 !-----------------------------------
2486 read (isccor,*,end=119,err=119) &
2487 nterm_sccor(i,j),nlor_sccor(i,j)
2491 do k=1,nterm_sccor(i,j)
2492 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2494 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2497 do k=1,nlor_sccor(i,j)
2498 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2499 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2500 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2501 (1+vlor3sccor(k,i,j)**2)
2503 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2512 write (iout,'(/a/)') 'Torsional constants:'
2515 write (iout,*) 'ityp',i,' jtyp',j
2516 write (iout,*) 'Fourier constants'
2517 do k=1,nterm_sccor(i,j)
2518 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2520 write (iout,*) 'Lorenz constants'
2521 do k=1,nlor_sccor(i,j)
2522 write (iout,'(3(1pe15.5))') &
2523 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2530 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2531 ! interaction energy of the Gly, Ala, and Pro prototypes.
2534 ! Read electrostatic-interaction parameters
2539 write (iout,'(/a)') 'Electrostatic interaction constants:'
2540 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2541 'IT','JT','APP','BPP','AEL6','AEL3'
2543 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2544 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2545 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2546 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2551 app (i,j)=epp(i,j)*rri*rri
2552 bpp (i,j)=-2.0D0*epp(i,j)*rri
2553 ael6(i,j)=elpp6(i,j)*4.2D0**6
2554 ael3(i,j)=elpp3(i,j)*4.2D0**3
2556 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2562 ! Read side-chain interaction parameters.
2564 !el from module energy - COMMON.INTERACT-------
2565 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2566 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2567 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2569 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2570 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2571 allocate(epslip(ntyp,ntyp))
2579 !--------------------------------
2581 read (isidep,*,end=117,err=117) ipot,expon
2582 if (ipot.lt.1 .or. ipot.gt.5) then
2583 write (iout,'(2a)') 'Error while reading SC interaction',&
2584 'potential file - unknown potential type.'
2586 call MPI_Finalize(Ierror)
2592 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2593 ', exponents are ',expon,2*expon
2594 ! goto (10,20,30,30,40) ipot
2596 !----------------------- LJ potential ---------------------------------
2598 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2599 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2600 (sigma0(i),i=1,ntyp)
2602 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2603 write (iout,'(a/)') 'The epsilon array:'
2604 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2605 write (iout,'(/a)') 'One-body parameters:'
2606 write (iout,'(a,4x,a)') 'residue','sigma'
2607 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2610 !----------------------- LJK potential --------------------------------
2612 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2613 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2614 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2616 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2617 write (iout,'(a/)') 'The epsilon array:'
2618 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2619 write (iout,'(/a)') 'One-body parameters:'
2620 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2621 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2625 !---------------------- GB or BP potential -----------------------------
2628 ! print *,"I AM in SCELE",scelemode
2629 if (scelemode.eq.0) then
2631 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2633 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2634 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2635 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2636 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2638 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2641 ! For the GB potential convert sigma'**2 into chi'
2644 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2648 write (iout,'(/a/)') 'Parameters of the BP potential:'
2649 write (iout,'(a/)') 'The epsilon array:'
2650 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2651 write (iout,'(/a)') 'One-body parameters:'
2652 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2654 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2655 chip(i),alp(i),i=1,ntyp)
2658 ! print *,ntyp,"NTYP"
2659 allocate(icharge(ntyp1))
2660 ! print *,ntyp,icharge(i)
2662 read (isidep,*) (icharge(i),i=1,ntyp)
2663 print *,ntyp,icharge(i)
2664 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2665 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2666 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2667 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2668 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2669 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2670 allocate(epsintab(ntyp,ntyp))
2671 allocate(dtail(2,ntyp,ntyp))
2672 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2673 allocate(wqdip(2,ntyp,ntyp))
2674 allocate(wstate(4,ntyp,ntyp))
2675 allocate(dhead(2,2,ntyp,ntyp))
2676 allocate(nstate(ntyp,ntyp))
2677 allocate(debaykap(ntyp,ntyp))
2679 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2680 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2684 ! write (*,*) "Im in ALAB", i, " ", j
2686 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2687 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2688 chis(i,j),chis(j,i), & !2 w tej linii
2689 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2690 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2691 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2692 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2693 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2694 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2695 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2696 IF ((LaTeX).and.(i.gt.24)) then
2697 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2698 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2699 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2700 chis(i,j),chis(j,i) !2 w tej linii
2702 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2707 IF ((LaTeX).and.(i.gt.24)) then
2708 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2709 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2710 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2711 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2712 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2713 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2714 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2721 sigma(i,j) = sigma(j,i)
2722 sigmap1(i,j)=sigmap1(j,i)
2723 sigmap2(i,j)=sigmap2(j,i)
2724 sigiso1(i,j)=sigiso1(j,i)
2725 sigiso2(i,j)=sigiso2(j,i)
2726 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2727 nstate(i,j) = nstate(j,i)
2728 dtail(1,i,j) = dtail(1,j,i)
2729 dtail(2,i,j) = dtail(2,j,i)
2731 alphasur(k,i,j) = alphasur(k,j,i)
2732 wstate(k,i,j) = wstate(k,j,i)
2733 alphiso(k,i,j) = alphiso(k,j,i)
2736 dhead(2,1,i,j) = dhead(1,1,j,i)
2737 dhead(2,2,i,j) = dhead(1,2,j,i)
2738 dhead(1,1,i,j) = dhead(2,1,j,i)
2739 dhead(1,2,i,j) = dhead(2,2,j,i)
2741 epshead(i,j) = epshead(j,i)
2742 sig0head(i,j) = sig0head(j,i)
2745 wqdip(k,i,j) = wqdip(k,j,i)
2748 wquad(i,j) = wquad(j,i)
2749 epsintab(i,j) = epsintab(j,i)
2750 debaykap(i,j)=debaykap(j,i)
2751 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2758 !--------------------- GBV potential -----------------------------------
2760 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2761 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2762 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2763 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2765 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2766 write (iout,'(a/)') 'The epsilon array:'
2767 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2768 write (iout,'(/a)') 'One-body parameters:'
2769 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2770 's||/s_|_^2',' chip ',' alph '
2771 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2772 sigii(i),chip(i),alp(i),i=1,ntyp)
2775 write(iout,*)"Wrong ipot"
2781 !-----------------------------------------------------------------------
2782 ! Calculate the "working" parameters of SC interactions.
2784 !el from module energy - COMMON.INTERACT-------
2785 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2786 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2787 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2788 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2789 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2790 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2792 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2798 if (scelemode.eq.0) then
2807 sc_aa_tube_par(:)=0.0d0
2808 sc_bb_tube_par(:)=0.0d0
2810 !--------------------------------
2815 epslip(i,j)=epslip(j,i)
2818 if (scelemode.eq.0) then
2821 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2822 sigma(j,i)=sigma(i,j)
2823 rs0(i,j)=dwa16*sigma(i,j)
2828 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2829 'Working parameters of the SC interactions:',&
2830 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2835 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2837 ! print *,"SIGMA ZLA?",sigma(i,j)
2845 sigeps=dsign(1.0D0,epsij)
2847 aa_aq(i,j)=epsij*rrij*rrij
2848 ! print *,"ADASKO",epsij,rrij,expon
2849 bb_aq(i,j)=-sigeps*epsij*rrij
2850 aa_aq(j,i)=aa_aq(i,j)
2851 bb_aq(j,i)=bb_aq(i,j)
2852 epsijlip=epslip(i,j)
2853 sigeps=dsign(1.0D0,epsijlip)
2854 epsijlip=dabs(epsijlip)
2855 aa_lip(i,j)=epsijlip*rrij*rrij
2856 bb_lip(i,j)=-sigeps*epsijlip*rrij
2857 aa_lip(j,i)=aa_lip(i,j)
2858 bb_lip(j,i)=bb_lip(i,j)
2859 !C write(iout,*) aa_lip
2860 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2861 sigt1sq=sigma0(i)**2
2862 sigt2sq=sigma0(j)**2
2865 ratsig1=sigt2sq/sigt1sq
2866 ratsig2=1.0D0/ratsig1
2867 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2868 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2869 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2873 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2874 sigmaii(i,j)=rsum_max
2875 sigmaii(j,i)=rsum_max
2877 ! sigmaii(i,j)=r0(i,j)
2878 ! sigmaii(j,i)=r0(i,j)
2880 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2881 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2882 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2883 augm(i,j)=epsij*r_augm**(2*expon)
2884 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2891 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2892 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2893 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2898 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2899 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2900 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2901 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2902 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2903 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2904 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2905 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2906 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2907 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2908 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2909 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2910 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2911 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2912 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2913 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2922 read (isidep_nucl,*) ipot_nucl
2923 ! print *,"TU?!",ipot_nucl
2924 if (ipot_nucl.eq.1) then
2925 do i=1,ntyp_molec(2)
2926 do j=i,ntyp_molec(2)
2927 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2928 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2932 do i=1,ntyp_molec(2)
2933 do j=i,ntyp_molec(2)
2934 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2935 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2936 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2940 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2941 do i=1,ntyp_molec(2)
2942 do j=i,ntyp_molec(2)
2943 rrij=sigma_nucl(i,j)
2947 epsij=4*eps_nucl(i,j)
2948 sigeps=dsign(1.0D0,epsij)
2950 aa_nucl(i,j)=epsij*rrij*rrij
2951 bb_nucl(i,j)=-sigeps*epsij*rrij
2952 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2953 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2954 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2955 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2956 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2957 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2960 aa_nucl(i,j)=aa_nucl(j,i)
2961 bb_nucl(i,j)=bb_nucl(j,i)
2962 ael3_nucl(i,j)=ael3_nucl(j,i)
2963 ael6_nucl(i,j)=ael6_nucl(j,i)
2964 ael63_nucl(i,j)=ael63_nucl(j,i)
2965 ael32_nucl(i,j)=ael32_nucl(j,i)
2966 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2967 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2968 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2969 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2970 eps_nucl(i,j)=eps_nucl(j,i)
2971 sigma_nucl(i,j)=sigma_nucl(j,i)
2972 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2976 write(iout,*) "tube param"
2977 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2978 ccavtubpep,dcavtubpep,tubetranenepep
2979 sigmapeptube=sigmapeptube**6
2980 sigeps=dsign(1.0D0,epspeptube)
2981 epspeptube=dabs(epspeptube)
2982 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2983 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2984 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2986 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2987 ccavtub(i),dcavtub(i),tubetranene(i)
2988 sigmasctube=sigmasctube**6
2989 sigeps=dsign(1.0D0,epssctube)
2990 epssctube=dabs(epssctube)
2991 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2992 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2993 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2995 !-----------------READING SC BASE POTENTIALS-----------------------------
2996 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2997 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2998 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2999 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3000 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3001 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3002 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3003 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3004 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3005 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3006 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3007 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3008 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3009 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3010 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3011 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3014 do i=1,ntyp_molec(1)
3015 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3016 write (*,*) "Im in ", i, " ", j
3017 read(isidep_scbase,*) &
3018 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3019 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3020 write(*,*) "eps",eps_scbase(i,j)
3021 read(isidep_scbase,*) &
3022 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3023 chis_scbase(i,j,1),chis_scbase(i,j,2)
3024 read(isidep_scbase,*) &
3025 dhead_scbasei(i,j), &
3026 dhead_scbasej(i,j), &
3027 rborn_scbasei(i,j),rborn_scbasej(i,j)
3028 read(isidep_scbase,*) &
3029 (wdipdip_scbase(k,i,j),k=1,3), &
3030 (wqdip_scbase(k,i,j),k=1,2)
3031 read(isidep_scbase,*) &
3032 alphapol_scbase(i,j), &
3033 epsintab_scbase(i,j)
3036 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3037 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3039 do i=1,ntyp_molec(1)
3040 do j=1,ntyp_molec(2)-1
3041 epsij=eps_scbase(i,j)
3042 rrij=sigma_scbase(i,j)
3047 sigeps=dsign(1.0D0,epsij)
3049 aa_scbase(i,j)=epsij*rrij*rrij
3050 bb_scbase(i,j)=-sigeps*epsij*rrij
3053 !-----------------READING PEP BASE POTENTIALS-------------------
3054 allocate(eps_pepbase(ntyp_molec(2)))
3055 allocate(sigma_pepbase(ntyp_molec(2)))
3056 allocate(chi_pepbase(ntyp_molec(2),2))
3057 allocate(chipp_pepbase(ntyp_molec(2),2))
3058 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3059 allocate(sigmap1_pepbase(ntyp_molec(2)))
3060 allocate(sigmap2_pepbase(ntyp_molec(2)))
3061 allocate(chis_pepbase(ntyp_molec(2),2))
3062 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3065 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3066 write (*,*) "Im in ", i, " ", j
3067 read(isidep_pepbase,*) &
3068 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3069 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3070 write(*,*) "eps",eps_pepbase(j)
3071 read(isidep_pepbase,*) &
3072 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3073 chis_pepbase(j,1),chis_pepbase(j,2)
3074 read(isidep_pepbase,*) &
3075 (wdipdip_pepbase(k,j),k=1,3)
3077 allocate(aa_pepbase(ntyp_molec(2)))
3078 allocate(bb_pepbase(ntyp_molec(2)))
3080 do j=1,ntyp_molec(2)-1
3081 epsij=eps_pepbase(j)
3082 rrij=sigma_pepbase(j)
3087 sigeps=dsign(1.0D0,epsij)
3089 aa_pepbase(j)=epsij*rrij*rrij
3090 bb_pepbase(j)=-sigeps*epsij*rrij
3092 !--------------READING SC PHOSPHATE-------------------------------------
3093 allocate(eps_scpho(ntyp_molec(1)))
3094 allocate(sigma_scpho(ntyp_molec(1)))
3095 allocate(chi_scpho(ntyp_molec(1),2))
3096 allocate(chipp_scpho(ntyp_molec(1),2))
3097 allocate(alphasur_scpho(4,ntyp_molec(1)))
3098 allocate(sigmap1_scpho(ntyp_molec(1)))
3099 allocate(sigmap2_scpho(ntyp_molec(1)))
3100 allocate(chis_scpho(ntyp_molec(1),2))
3101 allocate(wqq_scpho(ntyp_molec(1)))
3102 allocate(wqdip_scpho(2,ntyp_molec(1)))
3103 allocate(alphapol_scpho(ntyp_molec(1)))
3104 allocate(epsintab_scpho(ntyp_molec(1)))
3105 allocate(dhead_scphoi(ntyp_molec(1)))
3106 allocate(rborn_scphoi(ntyp_molec(1)))
3107 allocate(rborn_scphoj(ntyp_molec(1)))
3108 allocate(alphi_scpho(ntyp_molec(1)))
3112 do j=1,ntyp_molec(1) ! without U then we will take T for U
3113 write (*,*) "Im in scpho ", i, " ", j
3114 read(isidep_scpho,*) &
3115 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3116 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3117 write(*,*) "eps",eps_scpho(j)
3118 read(isidep_scpho,*) &
3119 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3120 chis_scpho(j,1),chis_scpho(j,2)
3121 read(isidep_scpho,*) &
3122 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3123 read(isidep_scpho,*) &
3124 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3128 allocate(aa_scpho(ntyp_molec(1)))
3129 allocate(bb_scpho(ntyp_molec(1)))
3131 do j=1,ntyp_molec(1)
3138 sigeps=dsign(1.0D0,epsij)
3140 aa_scpho(j)=epsij*rrij*rrij
3141 bb_scpho(j)=-sigeps*epsij*rrij
3145 read(isidep_peppho,*) &
3146 eps_peppho,sigma_peppho
3147 read(isidep_peppho,*) &
3148 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3149 read(isidep_peppho,*) &
3150 (wqdip_peppho(k),k=1,2)
3158 sigeps=dsign(1.0D0,epsij)
3160 aa_peppho=epsij*rrij*rrij
3161 bb_peppho=-sigeps*epsij*rrij
3164 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3169 ! Define the SC-p interaction constants (hard-coded; old style)
3172 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3174 ! aad(i,1)=0.3D0*4.0D0**12
3175 ! Following line for constants currently implemented
3176 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3177 aad(i,1)=1.5D0*4.0D0**12
3178 ! aad(i,1)=0.17D0*5.6D0**12
3180 ! "Soft" SC-p repulsion
3182 ! Following line for constants currently implemented
3183 ! aad(i,1)=0.3D0*4.0D0**6
3184 ! "Hard" SC-p repulsion
3185 bad(i,1)=3.0D0*4.0D0**6
3186 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3195 ! 8/9/01 Read the SC-p interaction constants from file
3198 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3201 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3202 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3203 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3204 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3208 write (iout,*) "Parameters of SC-p interactions:"
3210 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3211 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3216 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3218 do i=1,ntyp_molec(2)
3219 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3221 do i=1,ntyp_molec(2)
3222 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3223 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3225 r0pp=1.12246204830937298142*5.16158
3231 ! Define the constants of the disulfide bridge
3235 ! Old arbitrary potential - commented out.
3240 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3241 ! energy surface of diethyl disulfide.
3242 ! A. Liwo and U. Kozlowska, 11/24/03
3260 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3261 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3262 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3263 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3264 allocate(epsintabcat(ntyp,ntyp))
3265 allocate(dtailcat(2,ntyp,ntyp))
3266 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3267 allocate(wqdipcat(2,ntyp,ntyp))
3268 allocate(wstatecat(4,ntyp,ntyp))
3269 allocate(dheadcat(2,2,ntyp,ntyp))
3270 allocate(nstatecat(ntyp,ntyp))
3271 allocate(debaykapcat(ntyp,ntyp))
3273 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3274 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3275 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3276 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3277 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3280 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
3281 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3282 if (oldion.eq.0) then
3283 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3284 allocate(icharge(1:ntyp1))
3285 read(iion,*) (icharge(i),i=1,ntyp)
3290 do i=1,ntyp_molec(5)
3291 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3292 print *,msc(i,5),restok(i,5)
3296 do j=1,ntyp_molec(5)
3298 ! do j=1,ntyp_molec(5)
3299 ! write (*,*) "Im in ALAB", i, " ", j
3301 epscat(i,j),sigmacat(i,j), &
3302 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3303 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3305 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3306 ! chiscat(i,j),chiscat(j,i), &
3307 chis1cat(i,j),chis2cat(i,j), &
3309 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3310 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3311 dtailcat(1,i,j),dtailcat(2,i,j), &
3312 epsheadcat(i,j),sig0headcat(i,j), &
3314 ! rborncat(i,j),rborncat(j,i),&
3315 rborn1cat(i,j),rborn2cat(i,j),&
3316 (wqdipcat(k,i,j),k=1,2), &
3317 alphapolcat(i,j),alphapolcat(j,i), &
3318 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3320 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3321 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3323 ! write (iout,*) 'i= ', i, ' j= ', j
3324 ! write (iout,*) 'epsi0= ', epscat(1,j)
3325 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3326 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3327 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3328 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3329 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3330 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3331 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3332 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3333 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3334 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3335 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3336 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3337 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3338 ! write (iout,*) 'a1= ', rborncat(j,1)
3339 ! write (iout,*) 'a2= ', rborncat(1,j)
3340 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3341 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3342 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3343 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3347 ! If ((i.eq.1).and.(j.eq.27)) then
3348 ! write (iout,*) 'SEP'
3349 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3350 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3355 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3357 do j=1,ntyp_molec(5)
3361 sigeps=dsign(1.0D0,epsij)
3363 aa_aq_cat(i,j)=epsij*rrij*rrij
3364 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3369 do j=1,ntyp_molec(5)
3371 write (iout,*) 'i= ', i, ' j= ', j
3372 write (iout,*) 'epsi0= ', epscat(i,j)
3373 write (iout,*) 'sigma0= ', sigmacat(i,j)
3374 write (iout,*) 'chi1= ', chi1cat(i,j)
3375 write (iout,*) 'chi1= ', chi2cat(i,j)
3376 write (iout,*) 'chip1= ', chipp1cat(1,j)
3377 write (iout,*) 'chip2= ', chipp2cat(1,j)
3378 write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3379 write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3380 write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3381 write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3382 write (iout,*) 'sig1= ', sigmap1cat(1,j)
3383 write (iout,*) 'sig2= ', sigmap2cat(1,j)
3384 write (iout,*) 'chis1= ', chis1cat(1,j)
3385 write (iout,*) 'chis1= ', chis2cat(1,j)
3386 write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
3387 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
3388 write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3389 write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
3390 write (iout,*) 'a1= ', rborn1cat(i,j)
3391 write (iout,*) 'a2= ', rborn2cat(i,j)
3392 write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3393 write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3394 write (iout,*) 'alphapol2= ', alphapolcat(j,1)
3395 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3396 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3397 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
3400 If ((i.eq.1).and.(j.eq.27)) then
3401 write (iout,*) 'SEP'
3402 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3403 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3413 write (iout,'(/a)') "Disulfide bridge parameters:"
3414 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3415 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3416 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3417 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3420 if (shield_mode.gt.0) then
3421 pi=4.0D0*datan(1.0D0)
3422 !C VSolvSphere the volume of solving sphere
3424 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3425 !C there will be no distinction between proline peptide group and normal peptide
3426 !C group in case of shielding parameters
3427 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3428 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3429 write (iout,*) VSolvSphere,VSolvSphere_div
3430 !C long axis of side chain
3432 long_r_sidechain(i)=vbldsc0(1,i)
3433 ! if (scelemode.eq.0) then
3434 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3435 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3437 ! short_r_sidechain(i)=sigma(i,i)
3439 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3446 111 write (iout,*) "Error reading bending energy parameters."
3448 112 write (iout,*) "Error reading rotamer energy parameters."
3450 113 write (iout,*) "Error reading torsional energy parameters."
3452 114 write (iout,*) "Error reading double torsional energy parameters."
3454 115 write (iout,*) &
3455 "Error reading cumulant (multibody energy) parameters."
3457 116 write (iout,*) "Error reading electrostatic energy parameters."
3459 117 write (iout,*) "Error reading side chain interaction parameters."
3461 118 write (iout,*) "Error reading SCp interaction parameters."
3463 119 write (iout,*) "Error reading SCCOR parameters"
3465 121 write (iout,*) "Error in Czybyshev parameters"
3468 call MPI_Finalize(Ierror)
3472 end subroutine parmread
3474 !-----------------------------------------------------------------------------
3476 !-----------------------------------------------------------------------------
3477 subroutine printmat(ldim,m,n,iout,key,a)
3480 character(len=3),dimension(n) :: key
3481 real(kind=8),dimension(ldim,n) :: a
3483 integer :: i,j,k,m,iout,nlim
3487 write (iout,1000) (key(k),k=i,nlim)
3489 1000 format (/5x,8(6x,a3))
3490 1020 format (/80(1h-)/)
3492 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3495 1010 format (a3,2x,8(f9.4))
3497 end subroutine printmat
3498 !-----------------------------------------------------------------------------
3500 !-----------------------------------------------------------------------------
3502 ! Read the PDB file and convert the peptide geometry into virtual-chain
3505 use energy_data, only: itype,istype
3509 ! use control, only: rescode,sugarcode
3510 ! implicit real*8 (a-h,o-z)
3511 ! include 'DIMENSIONS'
3512 ! include 'COMMON.LOCAL'
3513 ! include 'COMMON.VAR'
3514 ! include 'COMMON.CHAIN'
3515 ! include 'COMMON.INTERACT'
3516 ! include 'COMMON.IOUNITS'
3517 ! include 'COMMON.GEO'
3518 ! include 'COMMON.NAMES'
3519 ! include 'COMMON.CONTROL'
3520 ! include 'COMMON.DISTFIT'
3521 ! include 'COMMON.SETUP'
3522 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3524 logical :: lprn=.true.,fail
3525 real(kind=8),dimension(3) :: e1,e2,e3
3526 real(kind=8) :: dcj,efree_temp
3527 character(len=3) :: seq,res,res2
3528 character(len=5) :: atom
3529 character(len=80) :: card
3530 real(kind=8),dimension(3,20) :: sccor
3531 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3532 integer :: isugar,molecprev,firstion
3533 character*1 :: sugar
3535 real(kind=8),dimension(3) :: ccc
3537 integer,dimension(2,maxres/3) :: hfrag_alloc
3538 integer,dimension(4,maxres/3) :: bfrag_alloc
3539 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3540 real(kind=8),dimension(:,:), allocatable :: c_temporary
3541 integer,dimension(:,:) , allocatable :: itype_temporary
3542 integer,dimension(:),allocatable :: istype_temp
3549 ! write (2,*) "UNRES_PDB",unres_pdb
3569 !-----------------------------
3570 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3571 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3572 if(.not. allocated(istype)) allocate(istype(maxres))
3574 read (ipdbin,'(a80)',end=10) card
3575 write (iout,'(a)') card
3576 if (card(:5).eq.'HELIX') then
3579 read(card(22:25),*) hfrag(1,nhfrag)
3580 read(card(34:37),*) hfrag(2,nhfrag)
3582 if (card(:5).eq.'SHEET') then
3585 read(card(24:26),*) bfrag(1,nbfrag)
3586 read(card(35:37),*) bfrag(2,nbfrag)
3587 !rc----------------------------------------
3588 !rc to be corrected !!!
3589 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3590 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3591 !rc----------------------------------------
3593 if (card(:3).eq.'END') then
3595 else if (card(:3).eq.'TER') then
3600 itype(ires_old,molecule)=ntyp1_molec(molecule)
3601 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3602 nres_molec(molecule)=nres_molec(molecule)+2
3604 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3607 dc(j,ires)=sccor(j,iii)
3610 call sccenter(ires,iii,sccor)
3616 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3617 ! Fish out the ATOM cards.
3618 ! write(iout,*) 'card',card(1:20)
3619 ! print *,"ATU ",card(1:6), CARD(3:6)
3621 if (index(card(1:4),'ATOM').gt.0) then
3622 read (card(12:16),*) atom
3623 ! write (iout,*) "! ",atom," !",ires
3624 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3625 read (card(23:26),*) ires
3626 read (card(18:20),'(a3)') res
3627 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3628 ! & " ires_old",ires_old
3629 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3630 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3631 if (ires-ishift+ishift1.ne.ires_old) then
3632 ! Calculate the CM of the preceding residue.
3633 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3635 ! write (iout,*) "Calculating sidechain center iii",iii
3638 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3641 call sccenter(ires_old,iii,sccor)
3645 ! Start new residue.
3646 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3649 else if (ibeg.eq.1) then
3650 write (iout,*) "BEG ires",ires
3652 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3655 nres_molec(molecule)=nres_molec(molecule)+1
3657 ires=ires-ishift+ishift1
3659 ! write (iout,*) "ishift",ishift," ires",ires,&
3660 ! " ires_old",ires_old
3662 else if (ibeg.eq.2) then
3664 ishift=-ires_old+ires-1 !!!!!
3665 ishift1=ishift1-1 !!!!!
3666 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3667 ires=ires-ishift+ishift1
3668 ! print *,ires,ishift,ishift1
3672 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3673 ires=ires-ishift+ishift1
3676 ! print *,'atom',ires,atom
3677 if (res.eq.'ACE' .or. res.eq.'NHE') then
3680 if (atom.eq.'CA '.or.atom.eq.'N ') then
3682 itype(ires,molecule)=rescode(ires,res,0,molecule)
3684 ! nres_molec(molecule)=nres_molec(molecule)+1
3688 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3689 ! nres_molec(molecule)=nres_molec(molecule)+1
3690 read (card(19:19),'(a1)') sugar
3691 isugar=sugarcode(sugar,ires)
3692 ! if (ibeg.eq.1) then
3696 ! print *,"ires=",ires,istype(ires)
3702 ires=ires-ishift+ishift1
3704 ! write (iout,*) "ires_old",ires_old," ires",ires
3705 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3708 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3709 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3710 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3711 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3712 ! print *,ires,ishift,ishift1
3713 ! write (iout,*) "backbone ",atom
3715 write (iout,'(2i3,2x,a,3f8.3)') &
3716 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3719 nres_molec(molecule)=nres_molec(molecule)+1
3721 sccor(j,iii)=c(j,ires)
3723 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3724 atom.eq."C2'" .or. atom.eq."C3'" &
3725 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3726 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3727 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3728 ! print *,ires,ishift,ishift1
3732 ! sccor(j,iii)=c(j,ires)
3735 c(j,ires)=c(j,ires)+ccc(j)/5.0
3737 print *,counter,molecule
3738 if (counter.eq.5) then
3740 nres_molec(molecule)=nres_molec(molecule)+1
3743 ! sccor(j,iii)=c(j,ires)
3747 ! print *, "ATOM",atom(1:3)
3748 else if (atom.eq."C5'") then
3749 read (card(19:19),'(a1)') sugar
3750 isugar=sugarcode(sugar,ires)
3755 ! print *,ires,istype(ires)
3759 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3760 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3761 nres_molec(molecule)=nres_molec(molecule)+1
3762 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3766 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3768 else if ((atom.eq."C1'").and.unres_pdb) then
3770 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3771 ! write (*,*) card(23:27),ires,itype(ires,1)
3772 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3773 atom.ne.'N' .and. atom.ne.'C' .and. &
3774 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3775 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3776 .and. atom.ne.'P '.and. &
3777 atom(1:1).ne.'H' .and. &
3778 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3780 ! write (iout,*) "sidechain ",atom
3781 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3782 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3783 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3785 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3788 ! print *,"IONS",ions,card(1:6)
3789 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3790 if (firstion.eq.0) then
3794 dc(j,ires)=sccor(j,iii)
3797 call sccenter(ires,iii,sccor)
3800 read (card(12:16),*) atom
3801 ! print *,"HETATOM", atom
3802 read (card(18:20),'(a3)') res
3803 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3804 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3805 .or.(atom(1:2).eq.'K ')) &
3808 if (molecule.ne.5) molecprev=molecule
3810 nres_molec(molecule)=nres_molec(molecule)+1
3811 print *,"HERE",nres_molec(molecule)
3813 itype(ires,molecule)=rescode(ires,res,0,molecule)
3814 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3818 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3819 if (ires.eq.0) return
3820 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3823 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3825 nres_molec(molecule)=nres_molec(molecule)-2
3826 print *,'I have',nres, nres_molec(:)
3828 do k=1,4 ! ions are without dummy
3829 if (nres_molec(k).eq.0) cycle
3831 ! write (iout,*) i,itype(i,1)
3832 ! if (itype(i,1).eq.ntyp1) then
3833 ! write (iout,*) "dummy",i,itype(i,1)
3835 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3836 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3840 if (itype(i,k).eq.ntyp1_molec(k)) then
3841 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3842 if (itype(i+2,k).eq.0) then
3843 ! print *,"masz sieczke"
3845 if (itype(i+2,j).ne.0) then
3847 itype(i+1,j)=ntyp1_molec(j)
3848 nres_molec(k)=nres_molec(k)-1
3849 nres_molec(j)=nres_molec(j)+1
3855 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3856 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3857 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3858 ! if (unres_pdb) then
3859 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3860 ! print *,i,'tu dochodze'
3861 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3869 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3873 dcj=(c(j,i-2)-c(j,i-3))/2.0
3874 if (dcj.eq.0) dcj=1.23591524223
3879 else !itype(i+1,1).eq.ntyp1
3880 ! if (unres_pdb) then
3881 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3882 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3889 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3890 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3894 dcj=(c(j,i+3)-c(j,i+2))/2.0
3895 if (dcj.eq.0) dcj=1.23591524223
3900 endif !itype(i+1,1).eq.ntyp1
3901 endif !itype.eq.ntyp1
3905 ! Calculate the CM of the last side chain.
3909 dc(j,ires)=sccor(j,iii)
3912 call sccenter(ires,iii,sccor)
3918 ! print *,"molecule",molecule
3919 if ((itype(nres,1).ne.10)) then
3921 if (molecule.eq.5) molecule=molecprev
3922 itype(nres,molecule)=ntyp1_molec(molecule)
3923 nres_molec(molecule)=nres_molec(molecule)+1
3924 ! if (unres_pdb) then
3925 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3926 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3933 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3937 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3938 c(j,nres)=c(j,nres-1)+dcj
3939 c(j,2*nres)=c(j,nres)
3943 ! print *,'I have',nres, nres_molec(:)
3945 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3947 if (nres.ne.nres0) then
3948 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3950 stop "Error nres value in WHAM input"
3953 !---------------------------------
3954 !el reallocate tables
3957 ! hfrag_alloc(j,i)=hfrag(j,i)
3960 ! bfrag_alloc(j,i)=bfrag(j,i)
3966 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3967 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3968 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3969 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3973 ! hfrag(j,i)=hfrag_alloc(j,i)
3978 ! bfrag(j,i)=bfrag_alloc(j,i)
3981 !el end reallocate tables
3982 !---------------------------------
3990 c(j,2*nres)=c(j,nres)
3993 if (itype(1,1).eq.ntyp1) then
3997 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3998 call refsys(2,3,4,e1,e2,e3,fail)
4005 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4006 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4010 dcj=(c(j,4)-c(j,3))/2.0
4016 ! First lets assign correct dummy to correct type of chain
4018 if (itype(1,1).eq.ntyp1) then
4019 if (itype(2,1).eq.0) then
4021 if (itype(2,j).ne.0) then
4023 itype(1,j)=ntyp1_molec(j)
4024 nres_molec(1)=nres_molec(1)-1
4025 nres_molec(j)=nres_molec(j)+1
4032 print *,'I have',nres, nres_molec(:)
4034 ! Copy the coordinates to reference coordinates
4040 ! Calculate internal coordinates.
4042 write (iout,'(/a)') &
4043 "Cartesian coordinates of the reference structure"
4044 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4045 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4047 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4048 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4049 (c(j,ires+nres),j=1,3)
4052 ! znamy już nres wiec mozna alokowac tablice
4053 ! Calculate internal coordinates.
4054 if(me.eq.king.or..not.out1file)then
4055 write (iout,'(a)') &
4056 "Backbone and SC coordinates as read from the PDB"
4058 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4059 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4060 (c(j,nres+ires),j=1,3)
4063 ! NOW LETS ROCK! SORTING
4064 allocate(c_temporary(3,2*nres))
4065 allocate(itype_temporary(nres,5))
4066 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4067 if (.not.allocated(istype)) write(iout,*) &
4068 "SOMETHING WRONG WITH ISTYTPE"
4069 allocate(istype_temp(nres))
4070 itype_temporary(:,:)=0
4074 if (itype(i,k).ne.0) then
4076 c_temporary(j,seqalingbegin)=c(j,i)
4077 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4080 itype_temporary(seqalingbegin,k)=itype(i,k)
4081 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4082 istype_temp(seqalingbegin)=istype(i)
4083 molnum(seqalingbegin)=k
4084 seqalingbegin=seqalingbegin+1
4090 c(j,i)=c_temporary(j,i)
4095 itype(i,k)=itype_temporary(i,k)
4096 istype(i)=istype_temp(i)
4099 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4100 ! I have only ions now dummy atoms in the system
4102 itype(1,5)=itype(2,5)
4105 itype(i,5)=itype(i+1,5)
4109 nres_molec(1)=nres_molec(1)-1
4111 ! if (itype(1,1).eq.ntyp1) then
4114 ! if (unres_pdb) then
4115 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4116 ! call refsys(2,3,4,e1,e2,e3,fail)
4123 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4127 ! dcj=(c(j,4)-c(j,3))/2.0
4129 ! c(j,nres+1)=c(j,1)
4135 write (iout,'(/a)') &
4136 "Cartesian coordinates of the reference structure after sorting"
4137 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4138 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4140 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4141 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4142 (c(j,ires+nres),j=1,3)
4146 ! print *,seqalingbegin,nres
4147 if(.not.allocated(vbld)) then
4148 allocate(vbld(2*nres))
4153 if(.not.allocated(vbld_inv)) then
4154 allocate(vbld_inv(2*nres))
4160 if(.not.allocated(theta)) then
4161 allocate(theta(nres+2))
4165 if(.not.allocated(phi)) allocate(phi(nres+2))
4166 if(.not.allocated(alph)) allocate(alph(nres+2))
4167 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4168 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4169 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4170 if(.not.allocated(costtab)) allocate(costtab(nres))
4171 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4172 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4173 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4174 if(.not.allocated(xxref)) allocate(xxref(nres))
4175 if(.not.allocated(yyref)) allocate(yyref(nres))
4176 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4177 if(.not.allocated(dc_norm)) then
4178 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4179 allocate(dc_norm(3,0:2*nres+2))
4183 call int_from_cart(.true.,.false.)
4184 call sc_loc_geom(.false.)
4186 thetaref(i)=theta(i)
4196 dc(j,i)=c(j,i+1)-c(j,i)
4197 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4202 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4203 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4205 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4209 ! Copy the coordinates to reference coordinates
4210 ! Splits to single chain if occurs
4214 ! cref(j,i,cou)=c(j,i)
4218 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4219 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4220 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4221 !-----------------------------
4225 write (iout,*) "symetr", symetr
4228 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4230 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4233 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4238 cref(j,i,cou)=c(j,i)
4239 cref(j,i+nres,cou)=c(j,i+nres)
4241 chain_rep(j,lll,kkk)=c(j,i)
4242 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4246 write (iout,*) chain_length
4247 if (chain_length.eq.0) chain_length=nres
4249 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4250 chain_rep(j,chain_length+nres,symetr) &
4251 =chain_rep(j,chain_length+nres,1)
4254 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4256 ! do kkk=1,chain_length
4257 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4261 ! makes copy of chains
4262 write (iout,*) "symetr", symetr
4267 if (symetr.gt.1) then
4274 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4280 write (iout,*) i,icha
4281 do lll=1,chain_length
4283 if (cou.le.nres) then
4285 kupa=mod(lll,chain_length)
4286 iprzes=(kkk-1)*chain_length+lll
4287 if (kupa.eq.0) kupa=chain_length
4288 write (iout,*) "kupa", kupa
4289 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4290 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4297 !-koniec robienia kopii
4300 write (iout,*) "nowa struktura", nperm
4302 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4304 cref(3,i,kkk),cref(1,nres+i,kkk),&
4305 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4307 100 format (//' alpha-carbon coordinates ',&
4308 ' centroid coordinates'/ &
4309 ' ', 6X,'X',11X,'Y',11X,'Z', &
4310 10X,'X',11X,'Y',11X,'Z')
4311 110 format (a,'(',i5,')',6f12.5)
4317 bfrag(i,j)=bfrag(i,j)-ishift
4323 hfrag(i,j)=hfrag(i,j)-ishift
4329 end subroutine readpdb
4330 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4331 !-----------------------------------------------------------------------------
4333 !-----------------------------------------------------------------------------
4334 subroutine read_control
4348 use random, only: random_init
4349 ! implicit real*8 (a-h,o-z)
4350 ! include 'DIMENSIONS'
4352 use prng, only:prng_restart
4354 logical :: OKRandom!, prng_restart
4357 ! include 'COMMON.IOUNITS'
4358 ! include 'COMMON.TIME1'
4359 ! include 'COMMON.THREAD'
4360 ! include 'COMMON.SBRIDGE'
4361 ! include 'COMMON.CONTROL'
4362 ! include 'COMMON.MCM'
4363 ! include 'COMMON.MAP'
4364 ! include 'COMMON.HEADER'
4365 ! include 'COMMON.CSA'
4366 ! include 'COMMON.CHAIN'
4367 ! include 'COMMON.MUCA'
4368 ! include 'COMMON.MD'
4369 ! include 'COMMON.FFIELD'
4370 ! include 'COMMON.INTERACT'
4371 ! include 'COMMON.SETUP'
4372 !el integer :: KDIAG,ICORFL,IXDR
4373 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4374 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4375 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4376 ! character(len=80) :: ucase
4377 character(len=640) :: controlcard
4379 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4385 read (INP,'(a)') titel
4386 call card_concat(controlcard,.true.)
4387 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4388 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4389 call reada(controlcard,'SEED',seed,0.0D0)
4390 call random_init(seed)
4391 ! Set up the time limit (caution! The time must be input in minutes!)
4392 read_cart=index(controlcard,'READ_CART').gt.0
4393 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4394 call readi(controlcard,'SYM',symetr,1)
4395 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4396 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4397 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4398 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4399 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4400 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4401 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4402 call reada(controlcard,'DRMS',drms,0.1D0)
4403 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4404 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4405 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4406 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4407 write (iout,'(a,f10.1)')'DRMS = ',drms
4408 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4409 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4411 call readi(controlcard,'NZ_START',nz_start,0)
4412 call readi(controlcard,'NZ_END',nz_end,0)
4413 call readi(controlcard,'IZ_SC',iz_sc,0)
4414 timlim=60.0D0*timlim
4415 safety = 60.0d0*safety
4418 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4419 !C SHIELD keyword sets if the shielding effect of side-chains is used
4420 !C 0 denots no shielding is used all peptide are equally despite the
4421 !C solvent accesible area
4422 !C 1 the newly introduced function
4423 !C 2 reseved for further possible developement
4424 call readi(controlcard,'SHIELD',shield_mode,0)
4425 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4426 write(iout,*) "shield_mode",shield_mode
4427 call readi(controlcard,'TORMODE',tor_mode,0)
4428 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4429 write(iout,*) "torsional and valence angle mode",tor_mode
4431 !C Varibles set size of box
4432 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4433 protein=index(controlcard,"PROTEIN").gt.0
4434 ions=index(controlcard,"IONS").gt.0
4435 call readi(controlcard,'OLDION',oldion,1)
4436 nucleic=index(controlcard,"NUCLEIC").gt.0
4437 write (iout,*) "with_theta_constr ",with_theta_constr
4438 AFMlog=(index(controlcard,'AFM'))
4439 selfguide=(index(controlcard,'SELFGUIDE'))
4440 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4441 call readi(controlcard,'GENCONSTR',genconstr,0)
4442 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4443 call reada(controlcard,'BOXY',boxysize,100.0d0)
4444 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4445 call readi(controlcard,'TUBEMOD',tubemode,0)
4446 print *,"SCELE",scelemode
4447 call readi(controlcard,"SCELEMODE",scelemode,0)
4448 print *,"SCELE",scelemode
4450 ! elemode = 0 is orignal UNRES electrostatics
4451 ! elemode = 1 is "Momo" potentials in progress
4452 ! elemode = 2 is in development EVALD
4455 write (iout,*) TUBEmode,"TUBEMODE"
4456 if (TUBEmode.gt.0) then
4457 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4458 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4459 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4460 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4461 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4462 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4463 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4464 buftubebot=bordtubebot+tubebufthick
4465 buftubetop=bordtubetop-tubebufthick
4468 ! CUTOFFF ON ELECTROSTATICS
4469 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4470 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4471 write(iout,*) "R_CUT_ELE=",r_cut_ele
4472 ! Lipidic parameters
4473 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4474 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4475 if (lipthick.gt.0.0d0) then
4476 bordliptop=(boxzsize+lipthick)/2.0
4477 bordlipbot=bordliptop-lipthick
4478 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4479 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4480 buflipbot=bordlipbot+lipbufthick
4481 bufliptop=bordliptop-lipbufthick
4482 if ((lipbufthick*2.0d0).gt.lipthick) &
4483 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4484 endif !lipthick.gt.0
4485 write(iout,*) "bordliptop=",bordliptop
4486 write(iout,*) "bordlipbot=",bordlipbot
4487 write(iout,*) "bufliptop=",bufliptop
4488 write(iout,*) "buflipbot=",buflipbot
4489 write (iout,*) "SHIELD MODE",shield_mode
4491 !C-------------------------
4492 minim=(index(controlcard,'MINIMIZE').gt.0)
4493 dccart=(index(controlcard,'CART').gt.0)
4494 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4495 overlapsc=.not.overlapsc
4496 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4497 searchsc=.not.searchsc
4498 sideadd=(index(controlcard,'SIDEADD').gt.0)
4499 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4500 outpdb=(index(controlcard,'PDBOUT').gt.0)
4501 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4502 pdbref=(index(controlcard,'PDBREF').gt.0)
4503 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4504 indpdb=index(controlcard,'PDBSTART')
4505 extconf=(index(controlcard,'EXTCONF').gt.0)
4506 call readi(controlcard,'IPRINT',iprint,0)
4507 call readi(controlcard,'MAXGEN',maxgen,10000)
4508 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4509 call readi(controlcard,"KDIAG",kdiag,0)
4510 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4511 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4512 write (iout,*) "RESCALE_MODE",rescale_mode
4513 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4514 if (index(controlcard,'REGULAR').gt.0.0D0) then
4515 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4519 if (index(controlcard,'CHECKGRAD').gt.0) then
4521 if (index(controlcard,'CART').gt.0) then
4523 elseif (index(controlcard,'CARINT').gt.0) then
4528 elseif (index(controlcard,'THREAD').gt.0) then
4530 call readi(controlcard,'THREAD',nthread,0)
4531 if (nthread.gt.0) then
4532 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4535 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4536 stop 'Error termination in Read_Control.'
4538 else if (index(controlcard,'MCMA').gt.0) then
4540 else if (index(controlcard,'MCEE').gt.0) then
4542 else if (index(controlcard,'MULTCONF').gt.0) then
4544 else if (index(controlcard,'MAP').gt.0) then
4546 call readi(controlcard,'MAP',nmap,0)
4547 else if (index(controlcard,'CSA').gt.0) then
4549 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4551 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4554 !fcm else if (index(controlcard,'MCMF').gt.0) then
4556 else if (index(controlcard,'SOFTREG').gt.0) then
4558 else if (index(controlcard,'CHECK_BOND').gt.0) then
4560 else if (index(controlcard,'TEST').gt.0) then
4562 else if (index(controlcard,'MD').gt.0) then
4564 else if (index(controlcard,'RE ').gt.0) then
4568 lmuca=index(controlcard,'MUCA').gt.0
4569 call readi(controlcard,'MUCADYN',mucadyn,0)
4570 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4571 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4573 write (iout,*) 'MUCADYN=',mucadyn
4574 write (iout,*) 'MUCASMOOTH=',muca_smooth
4577 iscode=index(controlcard,'ONE_LETTER')
4578 indphi=index(controlcard,'PHI')
4579 indback=index(controlcard,'BACK')
4580 iranconf=index(controlcard,'RAND_CONF')
4581 i2ndstr=index(controlcard,'USE_SEC_PRED')
4582 gradout=index(controlcard,'GRADOUT').gt.0
4583 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4584 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4585 if (me.eq.king .or. .not.out1file ) &
4586 write (iout,*) "DISTCHAINMAX",distchainmax
4588 if(me.eq.king.or..not.out1file) &
4589 write (iout,'(2a)') diagmeth(kdiag),&
4590 ' routine used to diagonalize matrices.'
4591 if (shield_mode.gt.0) then
4592 pi=4.0D0*datan(1.0D0)
4593 !C VSolvSphere the volume of solving sphere
4595 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4596 !C there will be no distinction between proline peptide group and normal peptide
4597 !C group in case of shielding parameters
4598 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4599 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4600 write (iout,*) VSolvSphere,VSolvSphere_div
4601 !C long axis of side chain
4603 ! long_r_sidechain(i)=vbldsc0(1,i)
4604 ! short_r_sidechain(i)=sigma0(i)
4605 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4612 end subroutine read_control
4613 !-----------------------------------------------------------------------------
4614 subroutine read_REMDpar
4616 ! Read REMD settings
4623 use control_data, only:out1file
4624 ! implicit real*8 (a-h,o-z)
4625 ! include 'DIMENSIONS'
4626 ! include 'COMMON.IOUNITS'
4627 ! include 'COMMON.TIME1'
4628 ! include 'COMMON.MD'
4631 !el include 'COMMON.LANGEVIN'
4633 !el include 'COMMON.LANGEVIN.lang0'
4635 ! include 'COMMON.INTERACT'
4636 ! include 'COMMON.NAMES'
4637 ! include 'COMMON.GEO'
4638 ! include 'COMMON.REMD'
4639 ! include 'COMMON.CONTROL'
4640 ! include 'COMMON.SETUP'
4641 ! character(len=80) :: ucase
4642 character(len=320) :: controlcard
4643 character(len=3200) :: controlcard1
4644 integer :: iremd_m_total
4647 ! real(kind=8) :: var,ene
4649 if(me.eq.king.or..not.out1file) &
4650 write (iout,*) "REMD setup"
4652 call card_concat(controlcard,.true.)
4653 call readi(controlcard,"NREP",nrep,3)
4654 call readi(controlcard,"NSTEX",nstex,1000)
4655 call reada(controlcard,"RETMIN",retmin,10.0d0)
4656 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4657 mremdsync=(index(controlcard,'SYNC').gt.0)
4658 call readi(controlcard,"NSYN",i_sync_step,100)
4659 restart1file=(index(controlcard,'REST1FILE').gt.0)
4660 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4661 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4662 if(max_cache_traj_use.gt.max_cache_traj) &
4663 max_cache_traj_use=max_cache_traj
4664 if(me.eq.king.or..not.out1file) then
4665 !d if (traj1file) then
4666 !rc caching is in testing - NTWX is not ignored
4667 !d write (iout,*) "NTWX value is ignored"
4668 !d write (iout,*) " trajectory is stored to one file by master"
4669 !d write (iout,*) " before exchange at NSTEX intervals"
4671 write (iout,*) "NREP= ",nrep
4672 write (iout,*) "NSTEX= ",nstex
4673 write (iout,*) "SYNC= ",mremdsync
4674 write (iout,*) "NSYN= ",i_sync_step
4675 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4678 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4679 if (index(controlcard,'TLIST').gt.0) then
4681 call card_concat(controlcard1,.true.)
4682 read(controlcard1,*) (remd_t(i),i=1,nrep)
4683 if(me.eq.king.or..not.out1file) &
4684 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4687 if (index(controlcard,'MLIST').gt.0) then
4689 call card_concat(controlcard1,.true.)
4690 read(controlcard1,*) (remd_m(i),i=1,nrep)
4691 if(me.eq.king.or..not.out1file) then
4692 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4695 iremd_m_total=iremd_m_total+remd_m(i)
4697 write (iout,*) 'Total number of replicas ',iremd_m_total
4700 if(me.eq.king.or..not.out1file) &
4701 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4703 end subroutine read_REMDpar
4704 !-----------------------------------------------------------------------------
4705 subroutine read_MDpar
4709 use control_data, only: r_cut,rlamb,out1file
4711 use geometry_data, only: pi
4713 ! implicit real*8 (a-h,o-z)
4714 ! include 'DIMENSIONS'
4715 ! include 'COMMON.IOUNITS'
4716 ! include 'COMMON.TIME1'
4717 ! include 'COMMON.MD'
4720 !el include 'COMMON.LANGEVIN'
4722 !el include 'COMMON.LANGEVIN.lang0'
4724 ! include 'COMMON.INTERACT'
4725 ! include 'COMMON.NAMES'
4726 ! include 'COMMON.GEO'
4727 ! include 'COMMON.SETUP'
4728 ! include 'COMMON.CONTROL'
4729 ! include 'COMMON.SPLITELE'
4730 ! character(len=80) :: ucase
4731 character(len=320) :: controlcard
4736 call card_concat(controlcard,.true.)
4737 call readi(controlcard,"NSTEP",n_timestep,1000000)
4738 call readi(controlcard,"NTWE",ntwe,100)
4739 call readi(controlcard,"NTWX",ntwx,1000)
4740 call reada(controlcard,"DT",d_time,1.0d-1)
4741 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4742 call reada(controlcard,"DAMAX",damax,1.0d1)
4743 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4744 call readi(controlcard,"LANG",lang,0)
4745 RESPA = index(controlcard,"RESPA") .gt. 0
4746 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4747 ntime_split0=ntime_split
4748 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4749 ntime_split0=ntime_split
4750 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4751 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4752 rest = index(controlcard,"REST").gt.0
4753 tbf = index(controlcard,"TBF").gt.0
4754 usampl = index(controlcard,"USAMPL").gt.0
4755 mdpdb = index(controlcard,"MDPDB").gt.0
4756 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4757 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4758 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4759 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4760 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4761 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4762 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4763 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4764 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4765 large = index(controlcard,"LARGE").gt.0
4766 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4767 rattle = index(controlcard,"RATTLE").gt.0
4768 preminim=(index(controlcard,'PREMINIM').gt.0)
4769 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
4770 write (iout,*) "PREMINIM ",preminim
4771 dccart=(index(controlcard,'CART').gt.0)
4772 if (preminim) call read_minim
4773 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4779 if(me.eq.king.or..not.out1file) then
4781 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4783 write (iout,'(a)') "The units are:"
4784 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4785 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4786 " acceleration: angstrom/(48.9 fs)**2"
4787 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4789 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4790 write (iout,'(a60,f10.5,a)') &
4791 "Initial time step of numerical integration:",d_time,&
4793 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4795 write (iout,'(2a,i4,a)') &
4796 "A-MTS algorithm used; initial time step for fast-varying",&
4797 " short-range forces split into",ntime_split," steps."
4798 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4799 r_cut," lambda",rlamb
4801 write (iout,'(2a,f10.5)') &
4802 "Maximum acceleration threshold to reduce the time step",&
4803 "/increase split number:",damax
4804 write (iout,'(2a,f10.5)') &
4805 "Maximum predicted energy drift to reduce the timestep",&
4806 "/increase split number:",edriftmax
4807 write (iout,'(a60,f10.5)') &
4808 "Maximum velocity threshold to reduce velocities:",dvmax
4809 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4810 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4811 if (rattle) write (iout,'(a60)') &
4812 "Rattle algorithm used to constrain the virtual bonds"
4816 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4817 call reada(controlcard,"RWAT",rwat,1.4d0)
4818 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4819 surfarea=index(controlcard,"SURFAREA").gt.0
4820 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4821 if(me.eq.king.or..not.out1file)then
4822 write (iout,'(/a,$)') "Langevin dynamics calculation"
4824 write (iout,'(a/)') &
4825 " with direct integration of Langevin equations"
4826 else if (lang.eq.2) then
4827 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4828 else if (lang.eq.3) then
4829 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4830 else if (lang.eq.4) then
4831 write (iout,'(a/)') " in overdamped mode"
4833 write (iout,'(//a,i5)') &
4834 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4837 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4838 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4839 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4840 write (iout,'(a60,f10.5)') &
4841 "Scaling factor of the friction forces:",scal_fric
4842 if (surfarea) write (iout,'(2a,i10,a)') &
4843 "Friction coefficients will be scaled by solvent-accessible",&
4844 " surface area every",reset_fricmat," steps."
4846 ! Calculate friction coefficients and bounds of stochastic forces
4847 eta=6*pi*cPoise*etawat
4848 if(me.eq.king.or..not.out1file) &
4849 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4852 do j=1,5 !types of molecules
4853 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4854 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4856 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4857 do j=1,5 !types of molecules
4859 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4860 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4864 if(me.eq.king.or..not.out1file)then
4865 write (iout,'(/2a/)') &
4866 "Radii of site types and friction coefficients and std's of",&
4867 " stochastic forces of fully exposed sites"
4868 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4870 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4871 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4875 if(me.eq.king.or..not.out1file)then
4876 write (iout,'(a)') "Berendsen bath calculation"
4877 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4878 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4880 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4881 count_reset_moment," steps"
4883 write (iout,'(a,i10,a)') &
4884 "Velocities will be reset at random every",count_reset_vel,&
4888 if(me.eq.king.or..not.out1file) &
4889 write (iout,'(a31)') "Microcanonical mode calculation"
4891 if(me.eq.king.or..not.out1file)then
4892 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4894 write(iout,*) "MD running with constraints."
4895 write(iout,*) "Equilibration time ", eq_time, " mtus."
4896 write(iout,*) "Constraining ", nfrag," fragments."
4897 write(iout,*) "Length of each fragment, weight and q0:"
4899 write (iout,*) "Set of restraints #",iset
4901 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4902 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4904 write(iout,*) "constraints between ", npair, "fragments."
4905 write(iout,*) "constraint pairs, weights and q0:"
4907 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4908 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4910 write(iout,*) "angle constraints within ", nfrag_back,&
4911 "backbone fragments."
4912 write(iout,*) "fragment, weights:"
4914 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4915 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4916 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4919 iset=mod(kolor,nset)+1
4922 if(me.eq.king.or..not.out1file) &
4923 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4925 end subroutine read_MDpar
4926 !-----------------------------------------------------------------------------
4930 ! implicit real*8 (a-h,o-z)
4931 ! include 'DIMENSIONS'
4932 ! include 'COMMON.MAP'
4933 ! include 'COMMON.IOUNITS'
4934 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4935 character(len=80) :: mapcard !,ucase
4938 ! real(kind=8) :: var,ene
4941 read (inp,'(a)') mapcard
4942 mapcard=ucase(mapcard)
4943 if (index(mapcard,'PHI').gt.0) then
4945 else if (index(mapcard,'THE').gt.0) then
4947 else if (index(mapcard,'ALP').gt.0) then
4949 else if (index(mapcard,'OME').gt.0) then
4952 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4953 stop 'Error - illegal variable spec in MAP card.'
4955 call readi (mapcard,'RES1',res1(imap),0)
4956 call readi (mapcard,'RES2',res2(imap),0)
4957 if (res1(imap).eq.0) then
4958 res1(imap)=res2(imap)
4959 else if (res2(imap).eq.0) then
4960 res2(imap)=res1(imap)
4962 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4963 write (iout,'(a)') &
4964 'Error - illegal definition of variable group in MAP.'
4965 stop 'Error - illegal definition of variable group in MAP.'
4967 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4968 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4969 call readi(mapcard,'NSTEP',nstep(imap),0)
4970 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4971 write (iout,'(a)') &
4972 'Illegal boundary and/or step size specification in MAP.'
4973 stop 'Illegal boundary and/or step size specification in MAP.'
4977 end subroutine map_read
4978 !-----------------------------------------------------------------------------
4981 use control_data, only: vdisulf
4983 ! implicit real*8 (a-h,o-z)
4984 ! include 'DIMENSIONS'
4985 ! include 'COMMON.IOUNITS'
4986 ! include 'COMMON.GEO'
4987 ! include 'COMMON.CSA'
4988 ! include 'COMMON.BANK'
4989 ! include 'COMMON.CONTROL'
4990 ! character(len=80) :: ucase
4991 character(len=620) :: mcmcard
4993 ! integer :: ntf,ik,iw_pdb
4994 ! real(kind=8) :: var,ene
4996 call card_concat(mcmcard,.true.)
4998 call readi(mcmcard,'NCONF',nconf,50)
4999 call readi(mcmcard,'NADD',nadd,0)
5000 call readi(mcmcard,'JSTART',jstart,1)
5001 call readi(mcmcard,'JEND',jend,1)
5002 call readi(mcmcard,'NSTMAX',nstmax,500000)
5003 call readi(mcmcard,'N0',n0,1)
5004 call readi(mcmcard,'N1',n1,6)
5005 call readi(mcmcard,'N2',n2,4)
5006 call readi(mcmcard,'N3',n3,0)
5007 call readi(mcmcard,'N4',n4,0)
5008 call readi(mcmcard,'N5',n5,0)
5009 call readi(mcmcard,'N6',n6,10)
5010 call readi(mcmcard,'N7',n7,0)
5011 call readi(mcmcard,'N8',n8,0)
5012 call readi(mcmcard,'N9',n9,0)
5013 call readi(mcmcard,'N14',n14,0)
5014 call readi(mcmcard,'N15',n15,0)
5015 call readi(mcmcard,'N16',n16,0)
5016 call readi(mcmcard,'N17',n17,0)
5017 call readi(mcmcard,'N18',n18,0)
5019 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5021 call readi(mcmcard,'NDIFF',ndiff,2)
5022 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5023 call readi(mcmcard,'IS1',is1,1)
5024 call readi(mcmcard,'IS2',is2,8)
5025 call readi(mcmcard,'NRAN0',nran0,4)
5026 call readi(mcmcard,'NRAN1',nran1,2)
5027 call readi(mcmcard,'IRR',irr,1)
5028 call readi(mcmcard,'NSEED',nseed,20)
5029 call readi(mcmcard,'NTOTAL',ntotal,10000)
5030 call reada(mcmcard,'CUT1',cut1,2.0d0)
5031 call reada(mcmcard,'CUT2',cut2,5.0d0)
5032 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5033 call readi(mcmcard,'ICMAX',icmax,3)
5034 call readi(mcmcard,'IRESTART',irestart,0)
5035 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5038 call reada(mcmcard,'DELE',dele,20.0d0)
5039 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5040 call readi(mcmcard,'IREF',iref,0)
5041 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5042 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5043 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5044 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5045 write (iout,*) "NCONF_IN",nconf_in
5047 end subroutine csaread
5048 !-----------------------------------------------------------------------------
5052 use control_data, only: MaxMoveType
5055 ! implicit real*8 (a-h,o-z)
5056 ! include 'DIMENSIONS'
5057 ! include 'COMMON.MCM'
5058 ! include 'COMMON.MCE'
5059 ! include 'COMMON.IOUNITS'
5060 ! character(len=80) :: ucase
5061 character(len=320) :: mcmcard
5064 ! real(kind=8) :: var,ene
5066 call card_concat(mcmcard,.true.)
5067 call readi(mcmcard,'MAXACC',maxacc,100)
5068 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5069 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5070 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5071 call readi(mcmcard,'MAXREPM',maxrepm,200)
5072 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5073 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5074 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5075 call reada(mcmcard,'E_UP',e_up,5.0D0)
5076 call reada(mcmcard,'DELTE',delte,0.1D0)
5077 call readi(mcmcard,'NSWEEP',nsweep,5)
5078 call readi(mcmcard,'NSTEPH',nsteph,0)
5079 call readi(mcmcard,'NSTEPC',nstepc,0)
5080 call reada(mcmcard,'TMIN',tmin,298.0D0)
5081 call reada(mcmcard,'TMAX',tmax,298.0D0)
5082 call readi(mcmcard,'NWINDOW',nwindow,0)
5083 call readi(mcmcard,'PRINT_MC',print_mc,0)
5084 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5085 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5086 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5087 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5088 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5089 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5090 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5091 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5092 if (nwindow.gt.0) then
5093 allocate(winstart(nwindow)) !!el (maxres)
5094 allocate(winend(nwindow)) !!el
5095 allocate(winlen(nwindow)) !!el
5096 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5098 winlen(i)=winend(i)-winstart(i)+1
5101 if (tmax.lt.tmin) tmax=tmin
5102 if (tmax.eq.tmin) then
5106 if (nstepc.gt.0 .and. nsteph.gt.0) then
5107 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5108 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5110 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5111 ! Probabilities of different move types
5112 sumpro_type(0)=0.0D0
5113 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5114 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5115 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5116 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5117 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5118 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5119 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5121 print *,'i',i,' sumprotype',sumpro_type(i)
5122 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5123 print *,'i',i,' sumprotype',sumpro_type(i)
5126 end subroutine mcmread
5127 !-----------------------------------------------------------------------------
5128 subroutine read_minim
5131 ! implicit real*8 (a-h,o-z)
5132 ! include 'DIMENSIONS'
5133 ! include 'COMMON.MINIM'
5134 ! include 'COMMON.IOUNITS'
5135 ! character(len=80) :: ucase
5136 character(len=320) :: minimcard
5138 ! integer :: ntf,ik,iw_pdb
5139 ! real(kind=8) :: var,ene
5141 call card_concat(minimcard,.true.)
5142 call readi(minimcard,'MAXMIN',maxmin,2000)
5143 call readi(minimcard,'MAXFUN',maxfun,5000)
5144 call readi(minimcard,'MINMIN',minmin,maxmin)
5145 call readi(minimcard,'MINFUN',minfun,maxmin)
5146 call reada(minimcard,'TOLF',tolf,1.0D-2)
5147 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5148 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5149 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5150 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5151 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5152 'Options in energy minimization:'
5153 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5154 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5155 'MinMin:',MinMin,' MinFun:',MinFun,&
5156 ' TolF:',TolF,' RTolF:',RTolF
5158 end subroutine read_minim
5159 !-----------------------------------------------------------------------------
5160 subroutine openunits
5162 use MD_data, only: usampl
5165 use control_data, only:out1file
5166 use control, only: getenv_loc
5167 ! implicit real*8 (a-h,o-z)
5168 ! include 'DIMENSIONS'
5171 character(len=16) :: form,nodename
5172 integer :: nodelen,ierror,npos
5174 ! include 'COMMON.SETUP'
5175 ! include 'COMMON.IOUNITS'
5176 ! include 'COMMON.MD'
5177 ! include 'COMMON.CONTROL'
5178 integer :: lenpre,lenpot,lentmp !,ilen
5180 character(len=3) :: out1file_text !,ucase
5181 character(len=3) :: ll
5184 ! integer :: ntf,ik,iw_pdb
5185 ! real(kind=8) :: var,ene
5187 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5188 call getenv_loc("PREFIX",prefix)
5190 call getenv_loc("POT",pot)
5191 call getenv_loc("DIRTMP",tmpdir)
5192 call getenv_loc("CURDIR",curdir)
5193 call getenv_loc("OUT1FILE",out1file_text)
5194 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5195 out1file_text=ucase(out1file_text)
5196 if (out1file_text(1:1).eq."Y") then
5199 out1file=fg_rank.gt.0
5204 if (lentmp.gt.0) then
5205 write (*,'(80(1h!))')
5206 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5207 write (*,'(80(1h!))')
5208 write (*,*)"All output files will be on node /tmp directory."
5210 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5211 if (me.eq.king) then
5212 write (*,*) "The master node is ",nodename
5213 else if (fg_rank.eq.0) then
5214 write (*,*) "I am the CG slave node ",nodename
5216 write (*,*) "I am the FG slave node ",nodename
5219 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5220 lenpre = lentmp+lenpre+1
5222 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5223 ! Get the names and open the input files
5224 #if defined(WINIFL) || defined(WINPGI)
5225 open(1,file=pref_orig(:ilen(pref_orig))// &
5226 '.inp',status='old',readonly,shared)
5227 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5228 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5229 ! Get parameter filenames and open the parameter files.
5230 call getenv_loc('BONDPAR',bondname)
5231 open (ibond,file=bondname,status='old',readonly,shared)
5232 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5233 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5234 call getenv_loc('THETPAR',thetname)
5235 open (ithep,file=thetname,status='old',readonly,shared)
5236 call getenv_loc('ROTPAR',rotname)
5237 open (irotam,file=rotname,status='old',readonly,shared)
5238 call getenv_loc('TORPAR',torname)
5239 open (itorp,file=torname,status='old',readonly,shared)
5240 call getenv_loc('TORDPAR',tordname)
5241 open (itordp,file=tordname,status='old',readonly,shared)
5242 call getenv_loc('FOURIER',fouriername)
5243 open (ifourier,file=fouriername,status='old',readonly,shared)
5244 call getenv_loc('ELEPAR',elename)
5245 open (ielep,file=elename,status='old',readonly,shared)
5246 call getenv_loc('SIDEPAR',sidename)
5247 open (isidep,file=sidename,status='old',readonly,shared)
5249 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5250 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5251 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5252 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5253 call getenv_loc('TORPAR_NUCL',torname_nucl)
5254 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5255 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5256 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5257 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5258 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5261 #elif (defined CRAY) || (defined AIX)
5262 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5264 ! print *,"Processor",myrank," opened file 1"
5265 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5266 ! print *,"Processor",myrank," opened file 9"
5267 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5268 ! Get parameter filenames and open the parameter files.
5269 call getenv_loc('BONDPAR',bondname)
5270 open (ibond,file=bondname,status='old',action='read')
5271 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5272 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5274 ! print *,"Processor",myrank," opened file IBOND"
5275 call getenv_loc('THETPAR',thetname)
5276 open (ithep,file=thetname,status='old',action='read')
5277 ! print *,"Processor",myrank," opened file ITHEP"
5278 call getenv_loc('ROTPAR',rotname)
5279 open (irotam,file=rotname,status='old',action='read')
5280 ! print *,"Processor",myrank," opened file IROTAM"
5281 call getenv_loc('TORPAR',torname)
5282 open (itorp,file=torname,status='old',action='read')
5283 ! print *,"Processor",myrank," opened file ITORP"
5284 call getenv_loc('TORDPAR',tordname)
5285 open (itordp,file=tordname,status='old',action='read')
5286 ! print *,"Processor",myrank," opened file ITORDP"
5287 call getenv_loc('SCCORPAR',sccorname)
5288 open (isccor,file=sccorname,status='old',action='read')
5289 ! print *,"Processor",myrank," opened file ISCCOR"
5290 call getenv_loc('FOURIER',fouriername)
5291 open (ifourier,file=fouriername,status='old',action='read')
5292 ! print *,"Processor",myrank," opened file IFOURIER"
5293 call getenv_loc('ELEPAR',elename)
5294 open (ielep,file=elename,status='old',action='read')
5295 ! print *,"Processor",myrank," opened file IELEP"
5296 call getenv_loc('SIDEPAR',sidename)
5297 open (isidep,file=sidename,status='old',action='read')
5299 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5300 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5301 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5302 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5303 call getenv_loc('TORPAR_NUCL',torname_nucl)
5304 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5305 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5306 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5307 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5308 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5310 call getenv_loc('LIPTRANPAR',liptranname)
5311 open (iliptranpar,file=liptranname,status='old',action='read')
5312 call getenv_loc('TUBEPAR',tubename)
5313 open (itube,file=tubename,status='old',action='read')
5314 call getenv_loc('IONPAR',ionname)
5315 open (iion,file=ionname,status='old',action='read')
5317 ! print *,"Processor",myrank," opened file ISIDEP"
5318 ! print *,"Processor",myrank," opened parameter files"
5320 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5321 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5322 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5323 ! Get parameter filenames and open the parameter files.
5324 call getenv_loc('BONDPAR',bondname)
5325 open (ibond,file=bondname,status='old')
5326 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5327 open (ibond_nucl,file=bondname_nucl,status='old')
5329 call getenv_loc('THETPAR',thetname)
5330 open (ithep,file=thetname,status='old')
5331 call getenv_loc('ROTPAR',rotname)
5332 open (irotam,file=rotname,status='old')
5333 call getenv_loc('TORPAR',torname)
5334 open (itorp,file=torname,status='old')
5335 call getenv_loc('TORDPAR',tordname)
5336 open (itordp,file=tordname,status='old')
5337 call getenv_loc('SCCORPAR',sccorname)
5338 open (isccor,file=sccorname,status='old')
5339 call getenv_loc('FOURIER',fouriername)
5340 open (ifourier,file=fouriername,status='old')
5341 call getenv_loc('ELEPAR',elename)
5342 open (ielep,file=elename,status='old')
5343 call getenv_loc('SIDEPAR',sidename)
5344 open (isidep,file=sidename,status='old')
5346 open (ithep_nucl,file=thetname_nucl,status='old')
5347 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5348 open (irotam_nucl,file=rotname_nucl,status='old')
5349 call getenv_loc('TORPAR_NUCL',torname_nucl)
5350 open (itorp_nucl,file=torname_nucl,status='old')
5351 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5352 open (itordp_nucl,file=tordname_nucl,status='old')
5353 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5354 open (isidep_nucl,file=sidename_nucl,status='old')
5356 call getenv_loc('LIPTRANPAR',liptranname)
5357 open (iliptranpar,file=liptranname,status='old')
5358 call getenv_loc('TUBEPAR',tubename)
5359 open (itube,file=tubename,status='old')
5360 call getenv_loc('IONPAR',ionname)
5361 open (iion,file=ionname,status='old')
5362 call getenv_loc('IONPAR_NUCL',ionnuclname)
5363 open (iionnucl,file=ionnuclname,status='old')
5365 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5367 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5368 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5369 ! Get parameter filenames and open the parameter files.
5370 call getenv_loc('BONDPAR',bondname)
5371 open (ibond,file=bondname,status='old',action='read')
5372 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5373 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5374 call getenv_loc('THETPAR',thetname)
5375 open (ithep,file=thetname,status='old',action='read')
5376 call getenv_loc('ROTPAR',rotname)
5377 open (irotam,file=rotname,status='old',action='read')
5378 call getenv_loc('TORPAR',torname)
5379 open (itorp,file=torname,status='old',action='read')
5380 call getenv_loc('TORDPAR',tordname)
5381 open (itordp,file=tordname,status='old',action='read')
5382 call getenv_loc('SCCORPAR',sccorname)
5383 open (isccor,file=sccorname,status='old',action='read')
5385 call getenv_loc('THETPARPDB',thetname_pdb)
5386 print *,"thetname_pdb ",thetname_pdb
5387 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5388 print *,ithep_pdb," opened"
5390 call getenv_loc('FOURIER',fouriername)
5391 open (ifourier,file=fouriername,status='old',readonly)
5392 call getenv_loc('ELEPAR',elename)
5393 open (ielep,file=elename,status='old',readonly)
5394 call getenv_loc('SIDEPAR',sidename)
5395 open (isidep,file=sidename,status='old',readonly)
5397 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5398 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5399 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5400 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5401 call getenv_loc('TORPAR_NUCL',torname_nucl)
5402 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5403 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5404 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5405 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5406 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5407 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5408 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5409 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5410 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5411 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5412 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5413 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5414 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5417 call getenv_loc('LIPTRANPAR',liptranname)
5418 open (iliptranpar,file=liptranname,status='old',action='read')
5419 call getenv_loc('TUBEPAR',tubename)
5420 open (itube,file=tubename,status='old',action='read')
5421 call getenv_loc('IONPAR',ionname)
5422 open (iion,file=ionname,status='old',action='read')
5423 call getenv_loc('IONPAR_NUCL',ionnuclname)
5424 open (iionnucl,file=ionnuclname,status='old',action='read')
5427 call getenv_loc('ROTPARPDB',rotname_pdb)
5428 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5431 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5432 #if defined(WINIFL) || defined(WINPGI)
5433 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5434 #elif (defined CRAY) || (defined AIX)
5435 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5437 open (iscpp_nucl,file=scpname_nucl,status='old')
5439 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5444 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5445 ! Use -DOLDSCP to use hard-coded constants instead.
5447 call getenv_loc('SCPPAR',scpname)
5448 #if defined(WINIFL) || defined(WINPGI)
5449 open (iscpp,file=scpname,status='old',readonly,shared)
5450 #elif (defined CRAY) || (defined AIX)
5451 open (iscpp,file=scpname,status='old',action='read')
5453 open (iscpp,file=scpname,status='old')
5455 open (iscpp,file=scpname,status='old',action='read')
5458 call getenv_loc('PATTERN',patname)
5459 #if defined(WINIFL) || defined(WINPGI)
5460 open (icbase,file=patname,status='old',readonly,shared)
5461 #elif (defined CRAY) || (defined AIX)
5462 open (icbase,file=patname,status='old',action='read')
5464 open (icbase,file=patname,status='old')
5466 open (icbase,file=patname,status='old',action='read')
5469 ! Open output file only for CG processes
5470 ! print *,"Processor",myrank," fg_rank",fg_rank
5471 if (fg_rank.eq.0) then
5473 if (nodes.eq.1) then
5476 npos = dlog10(dfloat(nodes-1))+1
5478 if (npos.lt.3) npos=3
5479 write (liczba,'(i1)') npos
5480 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5482 write (liczba,form) me
5483 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5484 liczba(:ilen(liczba))
5485 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5487 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5489 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5490 liczba(:ilen(liczba))//'.mol2'
5491 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5492 liczba(:ilen(liczba))//'.stat'
5494 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5495 //liczba(:ilen(liczba))//'.stat')
5496 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5499 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5500 liczba(:ilen(liczba))//'.const'
5505 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5506 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5507 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5508 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5509 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5511 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5513 rest2name=prefix(:ilen(prefix))//'.rst'
5515 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5518 #if defined(AIX) || defined(PGI)
5519 if (me.eq.king .or. .not. out1file) &
5520 open(iout,file=outname,status='unknown')
5522 if (fg_rank.gt.0) then
5523 write (liczba,'(i3.3)') myrank/nfgtasks
5524 write (ll,'(bz,i3.3)') fg_rank
5525 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5530 open(igeom,file=intname,status='unknown',position='append')
5531 open(ipdb,file=pdbname,status='unknown')
5532 open(imol2,file=mol2name,status='unknown')
5533 open(istat,file=statname,status='unknown',position='append')
5535 !1out open(iout,file=outname,status='unknown')
5538 if (me.eq.king .or. .not.out1file) &
5539 open(iout,file=outname,status='unknown')
5541 if (fg_rank.gt.0) then
5542 write (liczba,'(i3.3)') myrank/nfgtasks
5543 write (ll,'(bz,i3.3)') fg_rank
5544 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5549 open(igeom,file=intname,status='unknown',access='append')
5550 open(ipdb,file=pdbname,status='unknown')
5551 open(imol2,file=mol2name,status='unknown')
5552 open(istat,file=statname,status='unknown',access='append')
5554 !1out open(iout,file=outname,status='unknown')
5557 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5558 csa_seed=prefix(:lenpre)//'.CSA.seed'
5559 csa_history=prefix(:lenpre)//'.CSA.history'
5560 csa_bank=prefix(:lenpre)//'.CSA.bank'
5561 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5562 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5563 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5564 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5565 csa_int=prefix(:lenpre)//'.int'
5566 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5567 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5568 csa_in=prefix(:lenpre)//'.CSA.in'
5569 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5572 write (iout,'(80(1h-))')
5573 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5574 write (iout,'(80(1h-))')
5575 write (iout,*) "Input file : ",&
5576 pref_orig(:ilen(pref_orig))//'.inp'
5577 write (iout,*) "Output file : ",&
5578 outname(:ilen(outname))
5580 write (iout,*) "Sidechain potential file : ",&
5581 sidename(:ilen(sidename))
5583 write (iout,*) "SCp potential file : ",&
5584 scpname(:ilen(scpname))
5586 write (iout,*) "Electrostatic potential file : ",&
5587 elename(:ilen(elename))
5588 write (iout,*) "Cumulant coefficient file : ",&
5589 fouriername(:ilen(fouriername))
5590 write (iout,*) "Torsional parameter file : ",&
5591 torname(:ilen(torname))
5592 write (iout,*) "Double torsional parameter file : ",&
5593 tordname(:ilen(tordname))
5594 write (iout,*) "SCCOR parameter file : ",&
5595 sccorname(:ilen(sccorname))
5596 write (iout,*) "Bond & inertia constant file : ",&
5597 bondname(:ilen(bondname))
5598 write (iout,*) "Bending parameter file : ",&
5599 thetname(:ilen(thetname))
5600 write (iout,*) "Rotamer parameter file : ",&
5601 rotname(:ilen(rotname))
5604 write (iout,*) "Thetpdb parameter file : ",&
5605 thetname_pdb(:ilen(thetname_pdb))
5608 write (iout,*) "Threading database : ",&
5609 patname(:ilen(patname))
5611 write (iout,*)" DIRTMP : ",&
5613 write (iout,'(80(1h-))')
5616 end subroutine openunits
5617 !-----------------------------------------------------------------------------
5620 use geometry_data, only: nres,dc
5622 ! implicit real*8 (a-h,o-z)
5623 ! include 'DIMENSIONS'
5624 ! include 'COMMON.CHAIN'
5625 ! include 'COMMON.IOUNITS'
5626 ! include 'COMMON.MD'
5629 ! real(kind=8) :: var,ene
5631 open(irest2,file=rest2name,status='unknown')
5632 read(irest2,*) totT,EK,potE,totE,t_bath
5635 ! AL 4/17/17: Now reading d_t(0,:) too
5637 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5640 ! AL 4/17/17: Now reading d_c(0,:) too
5642 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5645 read (irest2,*) iset
5649 end subroutine readrst
5650 !-----------------------------------------------------------------------------
5651 subroutine read_fragments
5655 use control_data, only:out1file
5658 ! implicit real*8 (a-h,o-z)
5659 ! include 'DIMENSIONS'
5663 ! include 'COMMON.SETUP'
5664 ! include 'COMMON.CHAIN'
5665 ! include 'COMMON.IOUNITS'
5666 ! include 'COMMON.MD'
5667 ! include 'COMMON.CONTROL'
5670 ! real(kind=8) :: var,ene
5672 read(inp,*) nset,nfrag,npair,nfrag_back
5674 !el from module energy
5675 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5676 if(.not.allocated(wfrag_back)) then
5677 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5678 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5680 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5681 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5683 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5684 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5687 if(me.eq.king.or..not.out1file) &
5688 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5689 " nfrag_back",nfrag_back
5691 read(inp,*) mset(iset)
5693 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5695 if(me.eq.king.or..not.out1file) &
5696 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5697 ifrag(2,i,iset), qinfrag(i,iset)
5700 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5702 if(me.eq.king.or..not.out1file) &
5703 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5704 ipair(2,i,iset), qinpair(i,iset)
5707 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5708 wfrag_back(3,i,iset),&
5709 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5710 if(me.eq.king.or..not.out1file) &
5711 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5712 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5716 end subroutine read_fragments
5717 !-----------------------------------------------------------------------------
5719 !-----------------------------------------------------------------------------
5723 ! implicit real*8 (a-h,o-z)
5724 ! include 'DIMENSIONS'
5725 ! include 'COMMON.CSA'
5726 ! include 'COMMON.BANK'
5727 ! include 'COMMON.IOUNITS'
5729 ! integer :: ntf,ik,iw_pdb
5730 ! real(kind=8) :: var,ene
5732 open(icsa_in,file=csa_in,status="old",err=100)
5733 read(icsa_in,*) nconf
5734 read(icsa_in,*) jstart,jend
5735 read(icsa_in,*) nstmax
5736 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5737 read(icsa_in,*) nran0,nran1,irr
5738 read(icsa_in,*) nseed
5739 read(icsa_in,*) ntotal,cut1,cut2
5740 read(icsa_in,*) estop
5741 read(icsa_in,*) icmax,irestart
5742 read(icsa_in,*) ntbankm,dele,difcut
5743 read(icsa_in,*) iref,rmscut,pnccut
5744 read(icsa_in,*) ndiff
5751 end subroutine csa_read
5752 !-----------------------------------------------------------------------------
5753 subroutine initial_write
5756 ! implicit real*8 (a-h,o-z)
5757 ! include 'DIMENSIONS'
5758 ! include 'COMMON.CSA'
5759 ! include 'COMMON.BANK'
5760 ! include 'COMMON.IOUNITS'
5762 ! integer :: ntf,ik,iw_pdb
5763 ! real(kind=8) :: var,ene
5765 open(icsa_seed,file=csa_seed,status="unknown")
5766 write(icsa_seed,*) "seed"
5768 #if defined(AIX) || defined(PGI)
5769 open(icsa_history,file=csa_history,status="unknown",&
5772 open(icsa_history,file=csa_history,status="unknown",&
5775 write(icsa_history,*) nconf
5776 write(icsa_history,*) jstart,jend
5777 write(icsa_history,*) nstmax
5778 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5779 write(icsa_history,*) nran0,nran1,irr
5780 write(icsa_history,*) nseed
5781 write(icsa_history,*) ntotal,cut1,cut2
5782 write(icsa_history,*) estop
5783 write(icsa_history,*) icmax,irestart
5784 write(icsa_history,*) ntbankm,dele,difcut
5785 write(icsa_history,*) iref,rmscut,pnccut
5786 write(icsa_history,*) ndiff
5788 write(icsa_history,*)
5791 open(icsa_bank1,file=csa_bank1,status="unknown")
5792 write(icsa_bank1,*) 0
5796 end subroutine initial_write
5797 !-----------------------------------------------------------------------------
5798 subroutine restart_write
5801 ! implicit real*8 (a-h,o-z)
5802 ! include 'DIMENSIONS'
5803 ! include 'COMMON.IOUNITS'
5804 ! include 'COMMON.CSA'
5805 ! include 'COMMON.BANK'
5807 ! integer :: ntf,ik,iw_pdb
5808 ! real(kind=8) :: var,ene
5810 #if defined(AIX) || defined(PGI)
5811 open(icsa_history,file=csa_history,position="append")
5813 open(icsa_history,file=csa_history,access="append")
5815 write(icsa_history,*)
5816 write(icsa_history,*) "This is restart"
5817 write(icsa_history,*)
5818 write(icsa_history,*) nconf
5819 write(icsa_history,*) jstart,jend
5820 write(icsa_history,*) nstmax
5821 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5822 write(icsa_history,*) nran0,nran1,irr
5823 write(icsa_history,*) nseed
5824 write(icsa_history,*) ntotal,cut1,cut2
5825 write(icsa_history,*) estop
5826 write(icsa_history,*) icmax,irestart
5827 write(icsa_history,*) ntbankm,dele,difcut
5828 write(icsa_history,*) iref,rmscut,pnccut
5829 write(icsa_history,*) ndiff
5830 write(icsa_history,*)
5831 write(icsa_history,*) "irestart is: ", irestart
5833 write(icsa_history,*)
5837 end subroutine restart_write
5838 !-----------------------------------------------------------------------------
5840 !-----------------------------------------------------------------------------
5841 subroutine write_pdb(npdb,titelloc,ee)
5843 ! implicit real*8 (a-h,o-z)
5844 ! include 'DIMENSIONS'
5845 ! include 'COMMON.IOUNITS'
5846 character(len=50) :: titelloc1
5847 character*(*) :: titelloc
5848 character(len=3) :: zahl
5849 character(len=5) :: liczba5
5851 integer :: npdb !,ilen
5855 ! real(kind=8) :: var,ene
5859 if (npdb.lt.1000) then
5860 call numstr(npdb,zahl)
5861 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5863 if (npdb.lt.10000) then
5864 write(liczba5,'(i1,i4)') 0,npdb
5866 write(liczba5,'(i5)') npdb
5868 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5870 call pdbout(ee,titelloc1,ipdb)
5873 end subroutine write_pdb
5874 !-----------------------------------------------------------------------------
5876 !-----------------------------------------------------------------------------
5877 subroutine write_thread_summary
5878 ! Thread the sequence through a database of known structures
5879 use control_data, only: refstr
5881 use energy_data, only: n_ene_comp
5883 ! implicit real*8 (a-h,o-z)
5884 ! include 'DIMENSIONS'
5886 use MPI_data !include 'COMMON.INFO'
5889 ! include 'COMMON.CONTROL'
5890 ! include 'COMMON.CHAIN'
5891 ! include 'COMMON.DBASE'
5892 ! include 'COMMON.INTERACT'
5893 ! include 'COMMON.VAR'
5894 ! include 'COMMON.THREAD'
5895 ! include 'COMMON.FFIELD'
5896 ! include 'COMMON.SBRIDGE'
5897 ! include 'COMMON.HEADER'
5898 ! include 'COMMON.NAMES'
5899 ! include 'COMMON.IOUNITS'
5900 ! include 'COMMON.TIME1'
5902 integer,dimension(maxthread) :: ip
5903 real(kind=8),dimension(0:n_ene) :: energia
5905 integer :: i,j,ii,jj,ipj,ik,kk,ist
5906 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5908 write (iout,'(30x,a/)') &
5909 ' *********** Summary threading statistics ************'
5910 write (iout,'(a)') 'Initial energies:'
5911 write (iout,'(a4,2x,a12,14a14,3a8)') &
5912 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5913 'RMSnat','NatCONT','NNCONT','RMS'
5914 ! Energy sort patterns
5919 enet=ener(n_ene-1,ip(i))
5922 if (ener(n_ene-1,ip(j)).lt.enet) then
5924 enet=ener(n_ene-1,ip(j))
5936 ist=nres_base(2,ii)+ipatt(2,i)
5938 energia(i)=ener0(kk,i)
5940 etot=ener0(n_ene_comp+1,i)
5941 rmsnat=ener0(n_ene_comp+2,i)
5942 rms=ener0(n_ene_comp+3,i)
5943 frac=ener0(n_ene_comp+4,i)
5944 frac_nn=ener0(n_ene_comp+5,i)
5947 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5948 i,str_nam(ii),ist+1,&
5949 (energia(print_order(kk)),kk=1,nprint_ene),&
5950 etot,rmsnat,frac,frac_nn,rms
5952 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5953 i,str_nam(ii),ist+1,&
5954 (energia(print_order(kk)),kk=1,nprint_ene),etot
5957 write (iout,'(//a)') 'Final energies:'
5958 write (iout,'(a4,2x,a12,17a14,3a8)') &
5959 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5960 'RMSnat','NatCONT','NNCONT','RMS'
5964 ist=nres_base(2,ii)+ipatt(2,i)
5966 energia(kk)=ener(kk,ik)
5968 etot=ener(n_ene_comp+1,i)
5969 rmsnat=ener(n_ene_comp+2,i)
5970 rms=ener(n_ene_comp+3,i)
5971 frac=ener(n_ene_comp+4,i)
5972 frac_nn=ener(n_ene_comp+5,i)
5973 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5974 i,str_nam(ii),ist+1,&
5975 (energia(print_order(kk)),kk=1,nprint_ene),&
5976 etot,rmsnat,frac,frac_nn,rms
5978 write (iout,'(/a/)') 'IEXAM array:'
5979 write (iout,'(i5)') nexcl
5981 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5983 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5984 'Max. time for threading step ',max_time_for_thread,&
5985 'Average time for threading step: ',ave_time_for_thread
5987 end subroutine write_thread_summary
5988 !-----------------------------------------------------------------------------
5989 subroutine write_stat_thread(ithread,ipattern,ist)
5991 use energy_data, only: n_ene_comp
5993 ! implicit real*8 (a-h,o-z)
5994 ! include "DIMENSIONS"
5995 ! include "COMMON.CONTROL"
5996 ! include "COMMON.IOUNITS"
5997 ! include "COMMON.THREAD"
5998 ! include "COMMON.FFIELD"
5999 ! include "COMMON.DBASE"
6000 ! include "COMMON.NAMES"
6001 real(kind=8),dimension(0:n_ene) :: energia
6003 integer :: ithread,ipattern,ist,i
6004 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6006 #if defined(AIX) || defined(PGI)
6007 open(istat,file=statname,position='append')
6009 open(istat,file=statname,access='append')
6012 energia(i)=ener(i,ithread)
6014 etot=ener(n_ene_comp+1,ithread)
6015 rmsnat=ener(n_ene_comp+2,ithread)
6016 rms=ener(n_ene_comp+3,ithread)
6017 frac=ener(n_ene_comp+4,ithread)
6018 frac_nn=ener(n_ene_comp+5,ithread)
6019 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6020 ithread,str_nam(ipattern),ist+1,&
6021 (energia(print_order(i)),i=1,nprint_ene),&
6022 etot,rmsnat,frac,frac_nn,rms
6025 end subroutine write_stat_thread
6026 !-----------------------------------------------------------------------------
6028 !-----------------------------------------------------------------------------
6029 end module io_config