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(alphapolcat2(ntyp,ntyp))
3262 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3263 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3264 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3265 allocate(epsintabcat(ntyp,ntyp))
3266 allocate(dtailcat(2,ntyp,ntyp))
3267 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3268 allocate(wqdipcat(2,ntyp,ntyp))
3269 allocate(wstatecat(4,ntyp,ntyp))
3270 allocate(dheadcat(2,2,ntyp,ntyp))
3271 allocate(nstatecat(ntyp,ntyp))
3272 allocate(debaykapcat(ntyp,ntyp))
3274 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3275 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3276 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3277 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3278 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3281 if (.not.allocated(ichargecat)) allocate (ichargecat(ntyp_molec(5)))
3282 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3283 if (oldion.eq.0) then
3284 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3285 allocate(icharge(1:ntyp1))
3286 read(iion,*) (icharge(i),i=1,ntyp)
3291 do i=1,ntyp_molec(5)
3292 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3293 print *,msc(i,5),restok(i,5)
3299 do j=1,ntyp_molec(5)
3301 ! do j=1,ntyp_molec(5)
3302 ! write (*,*) "Im in ALAB", i, " ", j
3304 epscat(i,j),sigmacat(i,j), &
3305 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3306 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3308 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3309 ! chiscat(i,j),chiscat(j,i), &
3310 chis1cat(i,j),chis2cat(i,j), &
3312 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3313 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3314 dtailcat(1,i,j),dtailcat(2,i,j), &
3315 epsheadcat(i,j),sig0headcat(i,j), &
3317 ! rborncat(i,j),rborncat(j,i),&
3318 rborn1cat(i,j),rborn2cat(i,j),&
3319 (wqdipcat(k,i,j),k=1,2), &
3320 alphapolcat(i,j),alphapolcat2(j,i), &
3321 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3323 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3324 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3326 ! write (iout,*) 'i= ', i, ' j= ', j
3327 ! write (iout,*) 'epsi0= ', epscat(1,j)
3328 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3329 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3330 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3331 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3332 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3333 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3334 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3335 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3336 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3337 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3338 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3339 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3340 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3341 ! write (iout,*) 'a1= ', rborncat(j,1)
3342 ! write (iout,*) 'a2= ', rborncat(1,j)
3343 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3344 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3345 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3346 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3350 ! If ((i.eq.1).and.(j.eq.27)) then
3351 ! write (iout,*) 'SEP'
3352 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3353 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3358 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3360 do j=1,ntyp_molec(5)
3364 sigeps=dsign(1.0D0,epsij)
3366 aa_aq_cat(i,j)=epsij*rrij*rrij
3367 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3372 do j=1,ntyp_molec(5)
3374 write (iout,*) 'i= ', i, ' j= ', j
3375 write (iout,*) 'epsi0= ', epscat(i,j)
3376 write (iout,*) 'sigma0= ', sigmacat(i,j)
3377 write (iout,*) 'chi1= ', chi1cat(i,j)
3378 write (iout,*) 'chi1= ', chi2cat(i,j)
3379 write (iout,*) 'chip1= ', chipp1cat(i,j)
3380 write (iout,*) 'chip2= ', chipp2cat(i,j)
3381 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3382 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3383 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3384 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3385 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3386 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3387 write (iout,*) 'chis1= ', chis1cat(i,j)
3388 write (iout,*) 'chis1= ', chis2cat(i,j)
3389 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3390 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3391 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3392 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3393 write (iout,*) 'a1= ', rborn1cat(i,j)
3394 write (iout,*) 'a2= ', rborn2cat(i,j)
3395 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3396 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3397 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3398 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3399 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3400 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3403 If ((i.eq.1).and.(j.eq.27)) then
3404 write (iout,*) 'SEP'
3405 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3406 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3416 write (iout,'(/a)') "Disulfide bridge parameters:"
3417 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3418 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3419 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3420 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3423 if (shield_mode.gt.0) then
3424 pi=4.0D0*datan(1.0D0)
3425 !C VSolvSphere the volume of solving sphere
3427 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3428 !C there will be no distinction between proline peptide group and normal peptide
3429 !C group in case of shielding parameters
3430 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3431 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3432 write (iout,*) VSolvSphere,VSolvSphere_div
3433 !C long axis of side chain
3435 long_r_sidechain(i)=vbldsc0(1,i)
3436 ! if (scelemode.eq.0) then
3437 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3438 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3440 ! short_r_sidechain(i)=sigma(i,i)
3442 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3449 111 write (iout,*) "Error reading bending energy parameters."
3451 112 write (iout,*) "Error reading rotamer energy parameters."
3453 113 write (iout,*) "Error reading torsional energy parameters."
3455 114 write (iout,*) "Error reading double torsional energy parameters."
3457 115 write (iout,*) &
3458 "Error reading cumulant (multibody energy) parameters."
3460 116 write (iout,*) "Error reading electrostatic energy parameters."
3462 117 write (iout,*) "Error reading side chain interaction parameters."
3464 118 write (iout,*) "Error reading SCp interaction parameters."
3466 119 write (iout,*) "Error reading SCCOR parameters"
3468 121 write (iout,*) "Error in Czybyshev parameters"
3471 call MPI_Finalize(Ierror)
3475 end subroutine parmread
3477 !-----------------------------------------------------------------------------
3479 !-----------------------------------------------------------------------------
3480 subroutine printmat(ldim,m,n,iout,key,a)
3483 character(len=3),dimension(n) :: key
3484 real(kind=8),dimension(ldim,n) :: a
3486 integer :: i,j,k,m,iout,nlim
3490 write (iout,1000) (key(k),k=i,nlim)
3492 1000 format (/5x,8(6x,a3))
3493 1020 format (/80(1h-)/)
3495 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3498 1010 format (a3,2x,8(f9.4))
3500 end subroutine printmat
3501 !-----------------------------------------------------------------------------
3503 !-----------------------------------------------------------------------------
3505 ! Read the PDB file and convert the peptide geometry into virtual-chain
3508 use energy_data, only: itype,istype
3512 ! use control, only: rescode,sugarcode
3513 ! implicit real*8 (a-h,o-z)
3514 ! include 'DIMENSIONS'
3515 ! include 'COMMON.LOCAL'
3516 ! include 'COMMON.VAR'
3517 ! include 'COMMON.CHAIN'
3518 ! include 'COMMON.INTERACT'
3519 ! include 'COMMON.IOUNITS'
3520 ! include 'COMMON.GEO'
3521 ! include 'COMMON.NAMES'
3522 ! include 'COMMON.CONTROL'
3523 ! include 'COMMON.DISTFIT'
3524 ! include 'COMMON.SETUP'
3525 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3527 logical :: lprn=.true.,fail
3528 real(kind=8),dimension(3) :: e1,e2,e3
3529 real(kind=8) :: dcj,efree_temp
3530 character(len=3) :: seq,res,res2
3531 character(len=5) :: atom
3532 character(len=80) :: card
3533 real(kind=8),dimension(3,20) :: sccor
3534 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3535 integer :: isugar,molecprev,firstion
3536 character*1 :: sugar
3538 real(kind=8),dimension(3) :: ccc
3540 integer,dimension(2,maxres/3) :: hfrag_alloc
3541 integer,dimension(4,maxres/3) :: bfrag_alloc
3542 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3543 real(kind=8),dimension(:,:), allocatable :: c_temporary
3544 integer,dimension(:,:) , allocatable :: itype_temporary
3545 integer,dimension(:),allocatable :: istype_temp
3552 ! write (2,*) "UNRES_PDB",unres_pdb
3572 !-----------------------------
3573 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3574 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3575 if(.not. allocated(istype)) allocate(istype(maxres))
3577 read (ipdbin,'(a80)',end=10) card
3578 write (iout,'(a)') card
3579 if (card(:5).eq.'HELIX') then
3582 read(card(22:25),*) hfrag(1,nhfrag)
3583 read(card(34:37),*) hfrag(2,nhfrag)
3585 if (card(:5).eq.'SHEET') then
3588 read(card(24:26),*) bfrag(1,nbfrag)
3589 read(card(35:37),*) bfrag(2,nbfrag)
3590 !rc----------------------------------------
3591 !rc to be corrected !!!
3592 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3593 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3594 !rc----------------------------------------
3596 if (card(:3).eq.'END') then
3598 else if (card(:3).eq.'TER') then
3603 itype(ires_old,molecule)=ntyp1_molec(molecule)
3604 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3605 nres_molec(molecule)=nres_molec(molecule)+2
3607 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3610 dc(j,ires)=sccor(j,iii)
3613 call sccenter(ires,iii,sccor)
3619 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3620 ! Fish out the ATOM cards.
3621 ! write(iout,*) 'card',card(1:20)
3622 ! print *,"ATU ",card(1:6), CARD(3:6)
3624 if (index(card(1:4),'ATOM').gt.0) then
3625 read (card(12:16),*) atom
3626 ! write (iout,*) "! ",atom," !",ires
3627 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3628 read (card(23:26),*) ires
3629 read (card(18:20),'(a3)') res
3630 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3631 ! & " ires_old",ires_old
3632 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3633 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3634 if (ires-ishift+ishift1.ne.ires_old) then
3635 ! Calculate the CM of the preceding residue.
3636 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3638 ! write (iout,*) "Calculating sidechain center iii",iii
3641 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3644 call sccenter(ires_old,iii,sccor)
3648 ! Start new residue.
3649 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3652 else if (ibeg.eq.1) then
3653 write (iout,*) "BEG ires",ires
3655 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3658 nres_molec(molecule)=nres_molec(molecule)+1
3660 ires=ires-ishift+ishift1
3662 ! write (iout,*) "ishift",ishift," ires",ires,&
3663 ! " ires_old",ires_old
3665 else if (ibeg.eq.2) then
3667 ishift=-ires_old+ires-1 !!!!!
3668 ishift1=ishift1-1 !!!!!
3669 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3670 ires=ires-ishift+ishift1
3671 ! print *,ires,ishift,ishift1
3675 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3676 ires=ires-ishift+ishift1
3679 ! print *,'atom',ires,atom
3680 if (res.eq.'ACE' .or. res.eq.'NHE') then
3683 if (atom.eq.'CA '.or.atom.eq.'N ') then
3685 itype(ires,molecule)=rescode(ires,res,0,molecule)
3687 ! nres_molec(molecule)=nres_molec(molecule)+1
3691 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3692 ! nres_molec(molecule)=nres_molec(molecule)+1
3693 read (card(19:19),'(a1)') sugar
3694 isugar=sugarcode(sugar,ires)
3695 ! if (ibeg.eq.1) then
3699 ! print *,"ires=",ires,istype(ires)
3705 ires=ires-ishift+ishift1
3707 ! write (iout,*) "ires_old",ires_old," ires",ires
3708 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3711 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3712 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3713 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3714 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3715 ! print *,ires,ishift,ishift1
3716 ! write (iout,*) "backbone ",atom
3718 write (iout,'(2i3,2x,a,3f8.3)') &
3719 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3722 nres_molec(molecule)=nres_molec(molecule)+1
3724 sccor(j,iii)=c(j,ires)
3726 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3727 atom.eq."C2'" .or. atom.eq."C3'" &
3728 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3729 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3730 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3731 ! print *,ires,ishift,ishift1
3735 ! sccor(j,iii)=c(j,ires)
3738 c(j,ires)=c(j,ires)+ccc(j)/5.0
3740 print *,counter,molecule
3741 if (counter.eq.5) then
3743 nres_molec(molecule)=nres_molec(molecule)+1
3746 ! sccor(j,iii)=c(j,ires)
3750 ! print *, "ATOM",atom(1:3)
3751 else if (atom.eq."C5'") then
3752 read (card(19:19),'(a1)') sugar
3753 isugar=sugarcode(sugar,ires)
3758 ! print *,ires,istype(ires)
3762 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3763 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3764 nres_molec(molecule)=nres_molec(molecule)+1
3765 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3769 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3771 else if ((atom.eq."C1'").and.unres_pdb) then
3773 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3774 ! write (*,*) card(23:27),ires,itype(ires,1)
3775 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3776 atom.ne.'N' .and. atom.ne.'C' .and. &
3777 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3778 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3779 .and. atom.ne.'P '.and. &
3780 atom(1:1).ne.'H' .and. &
3781 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3783 ! write (iout,*) "sidechain ",atom
3784 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3785 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3786 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3788 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3791 ! print *,"IONS",ions,card(1:6)
3792 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3793 if (firstion.eq.0) then
3797 dc(j,ires)=sccor(j,iii)
3800 call sccenter(ires,iii,sccor)
3803 read (card(12:16),*) atom
3804 ! print *,"HETATOM", atom
3805 read (card(18:20),'(a3)') res
3806 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3807 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3808 .or.(atom(1:2).eq.'K ')) &
3811 if (molecule.ne.5) molecprev=molecule
3813 nres_molec(molecule)=nres_molec(molecule)+1
3814 print *,"HERE",nres_molec(molecule)
3816 itype(ires,molecule)=rescode(ires,res,0,molecule)
3817 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3821 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3822 if (ires.eq.0) return
3823 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3826 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3828 nres_molec(molecule)=nres_molec(molecule)-2
3829 print *,'I have',nres, nres_molec(:)
3831 do k=1,4 ! ions are without dummy
3832 if (nres_molec(k).eq.0) cycle
3834 ! write (iout,*) i,itype(i,1)
3835 ! if (itype(i,1).eq.ntyp1) then
3836 ! write (iout,*) "dummy",i,itype(i,1)
3838 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3839 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3843 if (itype(i,k).eq.ntyp1_molec(k)) then
3844 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3845 if (itype(i+2,k).eq.0) then
3846 ! print *,"masz sieczke"
3848 if (itype(i+2,j).ne.0) then
3850 itype(i+1,j)=ntyp1_molec(j)
3851 nres_molec(k)=nres_molec(k)-1
3852 nres_molec(j)=nres_molec(j)+1
3858 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3859 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3860 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3861 ! if (unres_pdb) then
3862 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3863 ! print *,i,'tu dochodze'
3864 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3872 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3876 dcj=(c(j,i-2)-c(j,i-3))/2.0
3877 if (dcj.eq.0) dcj=1.23591524223
3882 else !itype(i+1,1).eq.ntyp1
3883 ! if (unres_pdb) then
3884 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3885 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3892 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3893 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3897 dcj=(c(j,i+3)-c(j,i+2))/2.0
3898 if (dcj.eq.0) dcj=1.23591524223
3903 endif !itype(i+1,1).eq.ntyp1
3904 endif !itype.eq.ntyp1
3908 ! Calculate the CM of the last side chain.
3912 dc(j,ires)=sccor(j,iii)
3915 call sccenter(ires,iii,sccor)
3921 ! print *,"molecule",molecule
3922 if ((itype(nres,1).ne.10)) then
3924 if (molecule.eq.5) molecule=molecprev
3925 itype(nres,molecule)=ntyp1_molec(molecule)
3926 nres_molec(molecule)=nres_molec(molecule)+1
3927 ! if (unres_pdb) then
3928 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3929 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3936 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3940 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
3941 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
3942 c(j,2*nres)=c(j,nres)
3946 ! print *,'I have',nres, nres_molec(:)
3948 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3950 if (nres.ne.nres0) then
3951 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3953 stop "Error nres value in WHAM input"
3956 !---------------------------------
3957 !el reallocate tables
3960 ! hfrag_alloc(j,i)=hfrag(j,i)
3963 ! bfrag_alloc(j,i)=bfrag(j,i)
3969 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3970 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3971 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3972 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3976 ! hfrag(j,i)=hfrag_alloc(j,i)
3981 ! bfrag(j,i)=bfrag_alloc(j,i)
3984 !el end reallocate tables
3985 !---------------------------------
3993 c(j,2*nres)=c(j,nres)
3996 if (itype(1,1).eq.ntyp1) then
4000 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4001 call refsys(2,3,4,e1,e2,e3,fail)
4008 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4009 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4013 dcj=(c(j,4)-c(j,3))/2.0
4019 ! First lets assign correct dummy to correct type of chain
4021 if (itype(1,1).eq.ntyp1) then
4022 if (itype(2,1).eq.0) then
4024 if (itype(2,j).ne.0) then
4026 itype(1,j)=ntyp1_molec(j)
4027 nres_molec(1)=nres_molec(1)-1
4028 nres_molec(j)=nres_molec(j)+1
4035 print *,'I have',nres, nres_molec(:)
4037 ! Copy the coordinates to reference coordinates
4043 ! Calculate internal coordinates.
4045 write (iout,'(/a)') &
4046 "Cartesian coordinates of the reference structure"
4047 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4048 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4050 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4051 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4052 (c(j,ires+nres),j=1,3)
4055 ! znamy już nres wiec mozna alokowac tablice
4056 ! Calculate internal coordinates.
4057 if(me.eq.king.or..not.out1file)then
4058 write (iout,'(a)') &
4059 "Backbone and SC coordinates as read from the PDB"
4061 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4062 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4063 (c(j,nres+ires),j=1,3)
4066 ! NOW LETS ROCK! SORTING
4067 allocate(c_temporary(3,2*nres))
4068 allocate(itype_temporary(nres,5))
4069 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4070 if (.not.allocated(istype)) write(iout,*) &
4071 "SOMETHING WRONG WITH ISTYTPE"
4072 allocate(istype_temp(nres))
4073 itype_temporary(:,:)=0
4077 if (itype(i,k).ne.0) then
4079 c_temporary(j,seqalingbegin)=c(j,i)
4080 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4083 itype_temporary(seqalingbegin,k)=itype(i,k)
4084 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4085 istype_temp(seqalingbegin)=istype(i)
4086 molnum(seqalingbegin)=k
4087 seqalingbegin=seqalingbegin+1
4093 c(j,i)=c_temporary(j,i)
4098 itype(i,k)=itype_temporary(i,k)
4099 istype(i)=istype_temp(i)
4102 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4103 ! I have only ions now dummy atoms in the system
4105 itype(1,5)=itype(2,5)
4108 itype(i,5)=itype(i+1,5)
4112 nres_molec(1)=nres_molec(1)-1
4114 ! if (itype(1,1).eq.ntyp1) then
4117 ! if (unres_pdb) then
4118 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4119 ! call refsys(2,3,4,e1,e2,e3,fail)
4126 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4130 ! dcj=(c(j,4)-c(j,3))/2.0
4132 ! c(j,nres+1)=c(j,1)
4138 write (iout,'(/a)') &
4139 "Cartesian coordinates of the reference structure after sorting"
4140 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4141 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4143 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4144 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4145 (c(j,ires+nres),j=1,3)
4149 print *,seqalingbegin,nres
4150 if(.not.allocated(vbld)) then
4151 allocate(vbld(2*nres))
4156 if(.not.allocated(vbld_inv)) then
4157 allocate(vbld_inv(2*nres))
4163 if(.not.allocated(theta)) then
4164 allocate(theta(nres+2))
4168 if(.not.allocated(phi)) allocate(phi(nres+2))
4169 if(.not.allocated(alph)) allocate(alph(nres+2))
4170 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4171 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4172 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4173 if(.not.allocated(costtab)) allocate(costtab(nres))
4174 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4175 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4176 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4177 if(.not.allocated(xxref)) allocate(xxref(nres))
4178 if(.not.allocated(yyref)) allocate(yyref(nres))
4179 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4180 if(.not.allocated(dc_norm)) then
4181 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4182 allocate(dc_norm(3,0:2*nres+2))
4185 write(iout,*) "before int_from_cart"
4186 call int_from_cart(.true.,.false.)
4187 call sc_loc_geom(.false.)
4188 write(iout,*) "after int_from_cart"
4192 thetaref(i)=theta(i)
4195 write(iout,*) "after thetaref"
4203 dc(j,i)=c(j,i+1)-c(j,i)
4204 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4209 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4210 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4212 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4216 ! Copy the coordinates to reference coordinates
4217 ! Splits to single chain if occurs
4221 ! cref(j,i,cou)=c(j,i)
4227 ! chomo(j,i,k)=c(j,i)
4230 ! write(iout,*) "after chomo"
4232 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4233 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4234 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4235 !-----------------------------
4239 write (iout,*) "symetr", symetr
4242 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4244 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4247 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4252 cref(j,i,cou)=c(j,i)
4253 cref(j,i+nres,cou)=c(j,i+nres)
4255 chain_rep(j,lll,kkk)=c(j,i)
4256 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4260 write (iout,*) chain_length
4261 if (chain_length.eq.0) chain_length=nres
4263 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4264 chain_rep(j,chain_length+nres,symetr) &
4265 =chain_rep(j,chain_length+nres,1)
4268 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4270 ! do kkk=1,chain_length
4271 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4275 ! makes copy of chains
4276 write (iout,*) "symetr", symetr
4281 if (symetr.gt.1) then
4288 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4294 write (iout,*) i,icha
4295 do lll=1,chain_length
4297 if (cou.le.nres) then
4299 kupa=mod(lll,chain_length)
4300 iprzes=(kkk-1)*chain_length+lll
4301 if (kupa.eq.0) kupa=chain_length
4302 write (iout,*) "kupa", kupa
4303 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4304 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4311 !-koniec robienia kopii
4314 write (iout,*) "nowa struktura", nperm
4316 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4318 cref(3,i,kkk),cref(1,nres+i,kkk),&
4319 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4321 100 format (//' alpha-carbon coordinates ',&
4322 ' centroid coordinates'/ &
4323 ' ', 6X,'X',11X,'Y',11X,'Z', &
4324 10X,'X',11X,'Y',11X,'Z')
4325 110 format (a,'(',i5,')',6f12.5)
4331 bfrag(i,j)=bfrag(i,j)-ishift
4337 hfrag(i,j)=hfrag(i,j)-ishift
4343 end subroutine readpdb
4344 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4345 !-----------------------------------------------------------------------------
4347 !-----------------------------------------------------------------------------
4348 subroutine read_control
4362 use random, only: random_init
4363 ! implicit real*8 (a-h,o-z)
4364 ! include 'DIMENSIONS'
4366 use prng, only:prng_restart
4368 logical :: OKRandom!, prng_restart
4371 ! include 'COMMON.IOUNITS'
4372 ! include 'COMMON.TIME1'
4373 ! include 'COMMON.THREAD'
4374 ! include 'COMMON.SBRIDGE'
4375 ! include 'COMMON.CONTROL'
4376 ! include 'COMMON.MCM'
4377 ! include 'COMMON.MAP'
4378 ! include 'COMMON.HEADER'
4379 ! include 'COMMON.CSA'
4380 ! include 'COMMON.CHAIN'
4381 ! include 'COMMON.MUCA'
4382 ! include 'COMMON.MD'
4383 ! include 'COMMON.FFIELD'
4384 ! include 'COMMON.INTERACT'
4385 ! include 'COMMON.SETUP'
4386 !el integer :: KDIAG,ICORFL,IXDR
4387 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4388 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4389 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4390 ! character(len=80) :: ucase
4391 character(len=640) :: controlcard
4393 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4395 usampl=.false. ! the default value of usample should be 0
4396 ! write(iout,*) "KURWA2",usampl
4400 read (INP,'(a)') titel
4401 call card_concat(controlcard,.true.)
4402 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4403 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4404 call reada(controlcard,'SEED',seed,0.0D0)
4405 call random_init(seed)
4406 ! Set up the time limit (caution! The time must be input in minutes!)
4407 read_cart=index(controlcard,'READ_CART').gt.0
4408 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4409 call readi(controlcard,'SYM',symetr,1)
4410 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4411 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4412 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4413 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4414 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4415 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4416 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4417 call reada(controlcard,'DRMS',drms,0.1D0)
4418 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4419 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4420 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4421 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4422 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4423 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4424 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4425 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4426 write (iout,'(a,f10.1)')'DRMS = ',drms
4427 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4428 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4430 call readi(controlcard,'NZ_START',nz_start,0)
4431 call readi(controlcard,'NZ_END',nz_end,0)
4432 call readi(controlcard,'IZ_SC',iz_sc,0)
4433 timlim=60.0D0*timlim
4434 safety = 60.0d0*safety
4437 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4438 !C SHIELD keyword sets if the shielding effect of side-chains is used
4439 !C 0 denots no shielding is used all peptide are equally despite the
4440 !C solvent accesible area
4441 !C 1 the newly introduced function
4442 !C 2 reseved for further possible developement
4443 call readi(controlcard,'SHIELD',shield_mode,0)
4444 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4445 write(iout,*) "shield_mode",shield_mode
4446 call readi(controlcard,'TORMODE',tor_mode,0)
4447 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4448 write(iout,*) "torsional and valence angle mode",tor_mode
4450 !C Varibles set size of box
4451 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4452 protein=index(controlcard,"PROTEIN").gt.0
4453 ions=index(controlcard,"IONS").gt.0
4454 call readi(controlcard,'OLDION',oldion,1)
4455 nucleic=index(controlcard,"NUCLEIC").gt.0
4456 write (iout,*) "with_theta_constr ",with_theta_constr
4457 AFMlog=(index(controlcard,'AFM'))
4458 selfguide=(index(controlcard,'SELFGUIDE'))
4459 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4460 call readi(controlcard,'GENCONSTR',genconstr,0)
4461 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4462 call reada(controlcard,'BOXY',boxysize,100.0d0)
4463 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4464 call readi(controlcard,'TUBEMOD',tubemode,0)
4465 print *,"SCELE",scelemode
4466 call readi(controlcard,"SCELEMODE",scelemode,0)
4467 print *,"SCELE",scelemode
4469 ! elemode = 0 is orignal UNRES electrostatics
4470 ! elemode = 1 is "Momo" potentials in progress
4471 ! elemode = 2 is in development EVALD
4474 write (iout,*) TUBEmode,"TUBEMODE"
4475 if (TUBEmode.gt.0) then
4476 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4477 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4478 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4479 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4480 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4481 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4482 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4483 buftubebot=bordtubebot+tubebufthick
4484 buftubetop=bordtubetop-tubebufthick
4487 ! CUTOFFF ON ELECTROSTATICS
4488 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4489 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4490 write(iout,*) "R_CUT_ELE=",r_cut_ele
4491 ! Lipidic parameters
4492 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4493 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4494 if (lipthick.gt.0.0d0) then
4495 bordliptop=(boxzsize+lipthick)/2.0
4496 bordlipbot=bordliptop-lipthick
4497 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4498 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4499 buflipbot=bordlipbot+lipbufthick
4500 bufliptop=bordliptop-lipbufthick
4501 if ((lipbufthick*2.0d0).gt.lipthick) &
4502 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4503 endif !lipthick.gt.0
4504 write(iout,*) "bordliptop=",bordliptop
4505 write(iout,*) "bordlipbot=",bordlipbot
4506 write(iout,*) "bufliptop=",bufliptop
4507 write(iout,*) "buflipbot=",buflipbot
4508 write (iout,*) "SHIELD MODE",shield_mode
4510 !C-------------------------
4511 minim=(index(controlcard,'MINIMIZE').gt.0)
4512 dccart=(index(controlcard,'CART').gt.0)
4513 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4514 overlapsc=.not.overlapsc
4515 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4516 searchsc=.not.searchsc
4517 sideadd=(index(controlcard,'SIDEADD').gt.0)
4518 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4519 outpdb=(index(controlcard,'PDBOUT').gt.0)
4520 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4521 pdbref=(index(controlcard,'PDBREF').gt.0)
4522 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4523 indpdb=index(controlcard,'PDBSTART')
4524 extconf=(index(controlcard,'EXTCONF').gt.0)
4525 call readi(controlcard,'IPRINT',iprint,0)
4526 call readi(controlcard,'MAXGEN',maxgen,10000)
4527 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4528 call readi(controlcard,"KDIAG",kdiag,0)
4529 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4530 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4531 write (iout,*) "RESCALE_MODE",rescale_mode
4532 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4533 if (index(controlcard,'REGULAR').gt.0.0D0) then
4534 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4538 if (index(controlcard,'CHECKGRAD').gt.0) then
4540 if (index(controlcard,'CART').gt.0) then
4542 elseif (index(controlcard,'CARINT').gt.0) then
4547 elseif (index(controlcard,'THREAD').gt.0) then
4549 call readi(controlcard,'THREAD',nthread,0)
4550 if (nthread.gt.0) then
4551 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4554 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4555 stop 'Error termination in Read_Control.'
4557 else if (index(controlcard,'MCMA').gt.0) then
4559 else if (index(controlcard,'MCEE').gt.0) then
4561 else if (index(controlcard,'MULTCONF').gt.0) then
4563 else if (index(controlcard,'MAP').gt.0) then
4565 call readi(controlcard,'MAP',nmap,0)
4566 else if (index(controlcard,'CSA').gt.0) then
4568 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4570 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4573 !fcm else if (index(controlcard,'MCMF').gt.0) then
4575 else if (index(controlcard,'SOFTREG').gt.0) then
4577 else if (index(controlcard,'CHECK_BOND').gt.0) then
4579 else if (index(controlcard,'TEST').gt.0) then
4581 else if (index(controlcard,'MD').gt.0) then
4583 else if (index(controlcard,'RE ').gt.0) then
4587 lmuca=index(controlcard,'MUCA').gt.0
4588 call readi(controlcard,'MUCADYN',mucadyn,0)
4589 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4590 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4592 write (iout,*) 'MUCADYN=',mucadyn
4593 write (iout,*) 'MUCASMOOTH=',muca_smooth
4596 iscode=index(controlcard,'ONE_LETTER')
4597 indphi=index(controlcard,'PHI')
4598 indback=index(controlcard,'BACK')
4599 iranconf=index(controlcard,'RAND_CONF')
4600 i2ndstr=index(controlcard,'USE_SEC_PRED')
4601 gradout=index(controlcard,'GRADOUT').gt.0
4602 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4603 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4604 if (me.eq.king .or. .not.out1file ) &
4605 write (iout,*) "DISTCHAINMAX",distchainmax
4607 if(me.eq.king.or..not.out1file) &
4608 write (iout,'(2a)') diagmeth(kdiag),&
4609 ' routine used to diagonalize matrices.'
4610 if (shield_mode.gt.0) then
4611 pi=4.0D0*datan(1.0D0)
4612 !C VSolvSphere the volume of solving sphere
4614 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4615 !C there will be no distinction between proline peptide group and normal peptide
4616 !C group in case of shielding parameters
4617 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4618 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4619 write (iout,*) VSolvSphere,VSolvSphere_div
4620 !C long axis of side chain
4622 ! long_r_sidechain(i)=vbldsc0(1,i)
4623 ! short_r_sidechain(i)=sigma0(i)
4624 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4631 end subroutine read_control
4632 !-----------------------------------------------------------------------------
4633 subroutine read_REMDpar
4635 ! Read REMD settings
4642 use control_data, only:out1file
4643 ! implicit real*8 (a-h,o-z)
4644 ! include 'DIMENSIONS'
4645 ! include 'COMMON.IOUNITS'
4646 ! include 'COMMON.TIME1'
4647 ! include 'COMMON.MD'
4650 !el include 'COMMON.LANGEVIN'
4652 !el include 'COMMON.LANGEVIN.lang0'
4654 ! include 'COMMON.INTERACT'
4655 ! include 'COMMON.NAMES'
4656 ! include 'COMMON.GEO'
4657 ! include 'COMMON.REMD'
4658 ! include 'COMMON.CONTROL'
4659 ! include 'COMMON.SETUP'
4660 ! character(len=80) :: ucase
4661 character(len=320) :: controlcard
4662 character(len=3200) :: controlcard1
4663 integer :: iremd_m_total
4666 ! real(kind=8) :: var,ene
4668 if(me.eq.king.or..not.out1file) &
4669 write (iout,*) "REMD setup"
4671 call card_concat(controlcard,.true.)
4672 call readi(controlcard,"NREP",nrep,3)
4673 call readi(controlcard,"NSTEX",nstex,1000)
4674 call reada(controlcard,"RETMIN",retmin,10.0d0)
4675 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4676 mremdsync=(index(controlcard,'SYNC').gt.0)
4677 call readi(controlcard,"NSYN",i_sync_step,100)
4678 restart1file=(index(controlcard,'REST1FILE').gt.0)
4679 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4680 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4681 if(max_cache_traj_use.gt.max_cache_traj) &
4682 max_cache_traj_use=max_cache_traj
4683 if(me.eq.king.or..not.out1file) then
4684 !d if (traj1file) then
4685 !rc caching is in testing - NTWX is not ignored
4686 !d write (iout,*) "NTWX value is ignored"
4687 !d write (iout,*) " trajectory is stored to one file by master"
4688 !d write (iout,*) " before exchange at NSTEX intervals"
4690 write (iout,*) "NREP= ",nrep
4691 write (iout,*) "NSTEX= ",nstex
4692 write (iout,*) "SYNC= ",mremdsync
4693 write (iout,*) "NSYN= ",i_sync_step
4694 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4697 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4698 if (index(controlcard,'TLIST').gt.0) then
4700 call card_concat(controlcard1,.true.)
4701 read(controlcard1,*) (remd_t(i),i=1,nrep)
4702 if(me.eq.king.or..not.out1file) &
4703 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4706 if (index(controlcard,'MLIST').gt.0) then
4708 call card_concat(controlcard1,.true.)
4709 read(controlcard1,*) (remd_m(i),i=1,nrep)
4710 if(me.eq.king.or..not.out1file) then
4711 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4714 iremd_m_total=iremd_m_total+remd_m(i)
4716 write (iout,*) 'Total number of replicas ',iremd_m_total
4719 if(me.eq.king.or..not.out1file) &
4720 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4722 end subroutine read_REMDpar
4723 !-----------------------------------------------------------------------------
4724 subroutine read_MDpar
4728 use control_data, only: r_cut,rlamb,out1file
4730 use geometry_data, only: pi
4732 ! implicit real*8 (a-h,o-z)
4733 ! include 'DIMENSIONS'
4734 ! include 'COMMON.IOUNITS'
4735 ! include 'COMMON.TIME1'
4736 ! include 'COMMON.MD'
4739 !el include 'COMMON.LANGEVIN'
4741 !el include 'COMMON.LANGEVIN.lang0'
4743 ! include 'COMMON.INTERACT'
4744 ! include 'COMMON.NAMES'
4745 ! include 'COMMON.GEO'
4746 ! include 'COMMON.SETUP'
4747 ! include 'COMMON.CONTROL'
4748 ! include 'COMMON.SPLITELE'
4749 ! character(len=80) :: ucase
4750 character(len=320) :: controlcard
4755 call card_concat(controlcard,.true.)
4756 call readi(controlcard,"NSTEP",n_timestep,1000000)
4757 call readi(controlcard,"NTWE",ntwe,100)
4758 call readi(controlcard,"NTWX",ntwx,1000)
4759 call reada(controlcard,"DT",d_time,1.0d-1)
4760 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4761 call reada(controlcard,"DAMAX",damax,1.0d1)
4762 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4763 call readi(controlcard,"LANG",lang,0)
4764 RESPA = index(controlcard,"RESPA") .gt. 0
4765 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4766 ntime_split0=ntime_split
4767 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4768 ntime_split0=ntime_split
4769 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4770 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4771 rest = index(controlcard,"REST").gt.0
4772 tbf = index(controlcard,"TBF").gt.0
4773 usampl = index(controlcard,"USAMPL").gt.0
4774 ! write(iout,*) "KURWA",usampl
4775 mdpdb = index(controlcard,"MDPDB").gt.0
4776 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4777 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4778 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4779 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4780 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4781 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4782 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4783 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4784 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4785 large = index(controlcard,"LARGE").gt.0
4786 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4787 rattle = index(controlcard,"RATTLE").gt.0
4788 preminim=(index(controlcard,'PREMINIM').gt.0)
4789 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
4790 write (iout,*) "PREMINIM ",preminim
4791 dccart=(index(controlcard,'CART').gt.0)
4792 if (preminim) call read_minim
4793 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4799 if(me.eq.king.or..not.out1file) then
4801 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4803 write (iout,'(a)') "The units are:"
4804 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4805 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4806 " acceleration: angstrom/(48.9 fs)**2"
4807 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4809 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4810 write (iout,'(a60,f10.5,a)') &
4811 "Initial time step of numerical integration:",d_time,&
4813 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4815 write (iout,'(2a,i4,a)') &
4816 "A-MTS algorithm used; initial time step for fast-varying",&
4817 " short-range forces split into",ntime_split," steps."
4818 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4819 r_cut," lambda",rlamb
4821 write (iout,'(2a,f10.5)') &
4822 "Maximum acceleration threshold to reduce the time step",&
4823 "/increase split number:",damax
4824 write (iout,'(2a,f10.5)') &
4825 "Maximum predicted energy drift to reduce the timestep",&
4826 "/increase split number:",edriftmax
4827 write (iout,'(a60,f10.5)') &
4828 "Maximum velocity threshold to reduce velocities:",dvmax
4829 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4830 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4831 if (rattle) write (iout,'(a60)') &
4832 "Rattle algorithm used to constrain the virtual bonds"
4836 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4837 call reada(controlcard,"RWAT",rwat,1.4d0)
4838 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4839 surfarea=index(controlcard,"SURFAREA").gt.0
4840 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4841 if(me.eq.king.or..not.out1file)then
4842 write (iout,'(/a,$)') "Langevin dynamics calculation"
4844 write (iout,'(a/)') &
4845 " with direct integration of Langevin equations"
4846 else if (lang.eq.2) then
4847 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4848 else if (lang.eq.3) then
4849 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4850 else if (lang.eq.4) then
4851 write (iout,'(a/)') " in overdamped mode"
4853 write (iout,'(//a,i5)') &
4854 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4857 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4858 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4859 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4860 write (iout,'(a60,f10.5)') &
4861 "Scaling factor of the friction forces:",scal_fric
4862 if (surfarea) write (iout,'(2a,i10,a)') &
4863 "Friction coefficients will be scaled by solvent-accessible",&
4864 " surface area every",reset_fricmat," steps."
4866 ! Calculate friction coefficients and bounds of stochastic forces
4867 eta=6*pi*cPoise*etawat
4868 if(me.eq.king.or..not.out1file) &
4869 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4872 do j=1,5 !types of molecules
4873 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4874 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4876 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4877 do j=1,5 !types of molecules
4879 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4880 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4884 if(me.eq.king.or..not.out1file)then
4885 write (iout,'(/2a/)') &
4886 "Radii of site types and friction coefficients and std's of",&
4887 " stochastic forces of fully exposed sites"
4888 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4890 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4891 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4895 if(me.eq.king.or..not.out1file)then
4896 write (iout,'(a)') "Berendsen bath calculation"
4897 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4898 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4900 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4901 count_reset_moment," steps"
4903 write (iout,'(a,i10,a)') &
4904 "Velocities will be reset at random every",count_reset_vel,&
4908 if(me.eq.king.or..not.out1file) &
4909 write (iout,'(a31)') "Microcanonical mode calculation"
4911 if(me.eq.king.or..not.out1file)then
4912 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4914 write(iout,*) "MD running with constraints."
4915 write(iout,*) "Equilibration time ", eq_time, " mtus."
4916 write(iout,*) "Constraining ", nfrag," fragments."
4917 write(iout,*) "Length of each fragment, weight and q0:"
4919 write (iout,*) "Set of restraints #",iset
4921 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4922 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4924 write(iout,*) "constraints between ", npair, "fragments."
4925 write(iout,*) "constraint pairs, weights and q0:"
4927 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4928 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4930 write(iout,*) "angle constraints within ", nfrag_back,&
4931 "backbone fragments."
4932 write(iout,*) "fragment, weights:"
4934 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4935 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4936 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4939 iset=mod(kolor,nset)+1
4942 if(me.eq.king.or..not.out1file) &
4943 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4945 end subroutine read_MDpar
4946 !-----------------------------------------------------------------------------
4950 ! implicit real*8 (a-h,o-z)
4951 ! include 'DIMENSIONS'
4952 ! include 'COMMON.MAP'
4953 ! include 'COMMON.IOUNITS'
4954 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4955 character(len=80) :: mapcard !,ucase
4958 ! real(kind=8) :: var,ene
4961 read (inp,'(a)') mapcard
4962 mapcard=ucase(mapcard)
4963 if (index(mapcard,'PHI').gt.0) then
4965 else if (index(mapcard,'THE').gt.0) then
4967 else if (index(mapcard,'ALP').gt.0) then
4969 else if (index(mapcard,'OME').gt.0) then
4972 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4973 stop 'Error - illegal variable spec in MAP card.'
4975 call readi (mapcard,'RES1',res1(imap),0)
4976 call readi (mapcard,'RES2',res2(imap),0)
4977 if (res1(imap).eq.0) then
4978 res1(imap)=res2(imap)
4979 else if (res2(imap).eq.0) then
4980 res2(imap)=res1(imap)
4982 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4983 write (iout,'(a)') &
4984 'Error - illegal definition of variable group in MAP.'
4985 stop 'Error - illegal definition of variable group in MAP.'
4987 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4988 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4989 call readi(mapcard,'NSTEP',nstep(imap),0)
4990 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4991 write (iout,'(a)') &
4992 'Illegal boundary and/or step size specification in MAP.'
4993 stop 'Illegal boundary and/or step size specification in MAP.'
4997 end subroutine map_read
4998 !-----------------------------------------------------------------------------
5001 use control_data, only: vdisulf
5003 ! implicit real*8 (a-h,o-z)
5004 ! include 'DIMENSIONS'
5005 ! include 'COMMON.IOUNITS'
5006 ! include 'COMMON.GEO'
5007 ! include 'COMMON.CSA'
5008 ! include 'COMMON.BANK'
5009 ! include 'COMMON.CONTROL'
5010 ! character(len=80) :: ucase
5011 character(len=620) :: mcmcard
5013 ! integer :: ntf,ik,iw_pdb
5014 ! real(kind=8) :: var,ene
5016 call card_concat(mcmcard,.true.)
5018 call readi(mcmcard,'NCONF',nconf,50)
5019 call readi(mcmcard,'NADD',nadd,0)
5020 call readi(mcmcard,'JSTART',jstart,1)
5021 call readi(mcmcard,'JEND',jend,1)
5022 call readi(mcmcard,'NSTMAX',nstmax,500000)
5023 call readi(mcmcard,'N0',n0,1)
5024 call readi(mcmcard,'N1',n1,6)
5025 call readi(mcmcard,'N2',n2,4)
5026 call readi(mcmcard,'N3',n3,0)
5027 call readi(mcmcard,'N4',n4,0)
5028 call readi(mcmcard,'N5',n5,0)
5029 call readi(mcmcard,'N6',n6,10)
5030 call readi(mcmcard,'N7',n7,0)
5031 call readi(mcmcard,'N8',n8,0)
5032 call readi(mcmcard,'N9',n9,0)
5033 call readi(mcmcard,'N14',n14,0)
5034 call readi(mcmcard,'N15',n15,0)
5035 call readi(mcmcard,'N16',n16,0)
5036 call readi(mcmcard,'N17',n17,0)
5037 call readi(mcmcard,'N18',n18,0)
5039 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5041 call readi(mcmcard,'NDIFF',ndiff,2)
5042 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5043 call readi(mcmcard,'IS1',is1,1)
5044 call readi(mcmcard,'IS2',is2,8)
5045 call readi(mcmcard,'NRAN0',nran0,4)
5046 call readi(mcmcard,'NRAN1',nran1,2)
5047 call readi(mcmcard,'IRR',irr,1)
5048 call readi(mcmcard,'NSEED',nseed,20)
5049 call readi(mcmcard,'NTOTAL',ntotal,10000)
5050 call reada(mcmcard,'CUT1',cut1,2.0d0)
5051 call reada(mcmcard,'CUT2',cut2,5.0d0)
5052 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5053 call readi(mcmcard,'ICMAX',icmax,3)
5054 call readi(mcmcard,'IRESTART',irestart,0)
5055 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5058 call reada(mcmcard,'DELE',dele,20.0d0)
5059 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5060 call readi(mcmcard,'IREF',iref,0)
5061 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5062 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5063 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5064 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5065 write (iout,*) "NCONF_IN",nconf_in
5067 end subroutine csaread
5068 !-----------------------------------------------------------------------------
5072 use control_data, only: MaxMoveType
5075 ! implicit real*8 (a-h,o-z)
5076 ! include 'DIMENSIONS'
5077 ! include 'COMMON.MCM'
5078 ! include 'COMMON.MCE'
5079 ! include 'COMMON.IOUNITS'
5080 ! character(len=80) :: ucase
5081 character(len=320) :: mcmcard
5084 ! real(kind=8) :: var,ene
5086 call card_concat(mcmcard,.true.)
5087 call readi(mcmcard,'MAXACC',maxacc,100)
5088 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5089 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5090 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5091 call readi(mcmcard,'MAXREPM',maxrepm,200)
5092 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5093 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5094 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5095 call reada(mcmcard,'E_UP',e_up,5.0D0)
5096 call reada(mcmcard,'DELTE',delte,0.1D0)
5097 call readi(mcmcard,'NSWEEP',nsweep,5)
5098 call readi(mcmcard,'NSTEPH',nsteph,0)
5099 call readi(mcmcard,'NSTEPC',nstepc,0)
5100 call reada(mcmcard,'TMIN',tmin,298.0D0)
5101 call reada(mcmcard,'TMAX',tmax,298.0D0)
5102 call readi(mcmcard,'NWINDOW',nwindow,0)
5103 call readi(mcmcard,'PRINT_MC',print_mc,0)
5104 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5105 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5106 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5107 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5108 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5109 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5110 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5111 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5112 if (nwindow.gt.0) then
5113 allocate(winstart(nwindow)) !!el (maxres)
5114 allocate(winend(nwindow)) !!el
5115 allocate(winlen(nwindow)) !!el
5116 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5118 winlen(i)=winend(i)-winstart(i)+1
5121 if (tmax.lt.tmin) tmax=tmin
5122 if (tmax.eq.tmin) then
5126 if (nstepc.gt.0 .and. nsteph.gt.0) then
5127 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5128 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5130 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5131 ! Probabilities of different move types
5132 sumpro_type(0)=0.0D0
5133 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5134 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5135 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5136 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5137 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5138 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5139 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5141 print *,'i',i,' sumprotype',sumpro_type(i)
5142 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5143 print *,'i',i,' sumprotype',sumpro_type(i)
5146 end subroutine mcmread
5147 !-----------------------------------------------------------------------------
5148 subroutine read_minim
5151 ! implicit real*8 (a-h,o-z)
5152 ! include 'DIMENSIONS'
5153 ! include 'COMMON.MINIM'
5154 ! include 'COMMON.IOUNITS'
5155 ! character(len=80) :: ucase
5156 character(len=320) :: minimcard
5158 ! integer :: ntf,ik,iw_pdb
5159 ! real(kind=8) :: var,ene
5161 call card_concat(minimcard,.true.)
5162 call readi(minimcard,'MAXMIN',maxmin,2000)
5163 call readi(minimcard,'MAXFUN',maxfun,5000)
5164 call readi(minimcard,'MINMIN',minmin,maxmin)
5165 call readi(minimcard,'MINFUN',minfun,maxmin)
5166 call reada(minimcard,'TOLF',tolf,1.0D-2)
5167 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5168 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5169 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5170 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5171 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5172 'Options in energy minimization:'
5173 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5174 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5175 'MinMin:',MinMin,' MinFun:',MinFun,&
5176 ' TolF:',TolF,' RTolF:',RTolF
5178 end subroutine read_minim
5179 !-----------------------------------------------------------------------------
5180 subroutine openunits
5182 use MD_data, only: usampl
5185 use control_data, only:out1file
5186 use control, only: getenv_loc
5187 ! implicit real*8 (a-h,o-z)
5188 ! include 'DIMENSIONS'
5191 character(len=16) :: form,nodename
5192 integer :: nodelen,ierror,npos
5194 ! include 'COMMON.SETUP'
5195 ! include 'COMMON.IOUNITS'
5196 ! include 'COMMON.MD'
5197 ! include 'COMMON.CONTROL'
5198 integer :: lenpre,lenpot,lentmp !,ilen
5200 character(len=3) :: out1file_text !,ucase
5201 character(len=3) :: ll
5204 ! integer :: ntf,ik,iw_pdb
5205 ! real(kind=8) :: var,ene
5207 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5208 call getenv_loc("PREFIX",prefix)
5210 call getenv_loc("POT",pot)
5211 call getenv_loc("DIRTMP",tmpdir)
5212 call getenv_loc("CURDIR",curdir)
5213 call getenv_loc("OUT1FILE",out1file_text)
5214 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5215 out1file_text=ucase(out1file_text)
5216 if (out1file_text(1:1).eq."Y") then
5219 out1file=fg_rank.gt.0
5224 if (lentmp.gt.0) then
5225 write (*,'(80(1h!))')
5226 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5227 write (*,'(80(1h!))')
5228 write (*,*)"All output files will be on node /tmp directory."
5230 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5231 if (me.eq.king) then
5232 write (*,*) "The master node is ",nodename
5233 else if (fg_rank.eq.0) then
5234 write (*,*) "I am the CG slave node ",nodename
5236 write (*,*) "I am the FG slave node ",nodename
5239 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5240 lenpre = lentmp+lenpre+1
5242 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5243 ! Get the names and open the input files
5244 #if defined(WINIFL) || defined(WINPGI)
5245 open(1,file=pref_orig(:ilen(pref_orig))// &
5246 '.inp',status='old',readonly,shared)
5247 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5248 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5249 ! Get parameter filenames and open the parameter files.
5250 call getenv_loc('BONDPAR',bondname)
5251 open (ibond,file=bondname,status='old',readonly,shared)
5252 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5253 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5254 call getenv_loc('THETPAR',thetname)
5255 open (ithep,file=thetname,status='old',readonly,shared)
5256 call getenv_loc('ROTPAR',rotname)
5257 open (irotam,file=rotname,status='old',readonly,shared)
5258 call getenv_loc('TORPAR',torname)
5259 open (itorp,file=torname,status='old',readonly,shared)
5260 call getenv_loc('TORDPAR',tordname)
5261 open (itordp,file=tordname,status='old',readonly,shared)
5262 call getenv_loc('FOURIER',fouriername)
5263 open (ifourier,file=fouriername,status='old',readonly,shared)
5264 call getenv_loc('ELEPAR',elename)
5265 open (ielep,file=elename,status='old',readonly,shared)
5266 call getenv_loc('SIDEPAR',sidename)
5267 open (isidep,file=sidename,status='old',readonly,shared)
5269 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5270 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5271 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5272 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5273 call getenv_loc('TORPAR_NUCL',torname_nucl)
5274 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5275 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5276 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5277 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5278 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5281 #elif (defined CRAY) || (defined AIX)
5282 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5284 ! print *,"Processor",myrank," opened file 1"
5285 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5286 ! print *,"Processor",myrank," opened file 9"
5287 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5288 ! Get parameter filenames and open the parameter files.
5289 call getenv_loc('BONDPAR',bondname)
5290 open (ibond,file=bondname,status='old',action='read')
5291 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5292 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5294 ! print *,"Processor",myrank," opened file IBOND"
5295 call getenv_loc('THETPAR',thetname)
5296 open (ithep,file=thetname,status='old',action='read')
5297 ! print *,"Processor",myrank," opened file ITHEP"
5298 call getenv_loc('ROTPAR',rotname)
5299 open (irotam,file=rotname,status='old',action='read')
5300 ! print *,"Processor",myrank," opened file IROTAM"
5301 call getenv_loc('TORPAR',torname)
5302 open (itorp,file=torname,status='old',action='read')
5303 ! print *,"Processor",myrank," opened file ITORP"
5304 call getenv_loc('TORDPAR',tordname)
5305 open (itordp,file=tordname,status='old',action='read')
5306 ! print *,"Processor",myrank," opened file ITORDP"
5307 call getenv_loc('SCCORPAR',sccorname)
5308 open (isccor,file=sccorname,status='old',action='read')
5309 ! print *,"Processor",myrank," opened file ISCCOR"
5310 call getenv_loc('FOURIER',fouriername)
5311 open (ifourier,file=fouriername,status='old',action='read')
5312 ! print *,"Processor",myrank," opened file IFOURIER"
5313 call getenv_loc('ELEPAR',elename)
5314 open (ielep,file=elename,status='old',action='read')
5315 ! print *,"Processor",myrank," opened file IELEP"
5316 call getenv_loc('SIDEPAR',sidename)
5317 open (isidep,file=sidename,status='old',action='read')
5319 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5320 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5321 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5322 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5323 call getenv_loc('TORPAR_NUCL',torname_nucl)
5324 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5325 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5326 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5327 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5328 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5330 call getenv_loc('LIPTRANPAR',liptranname)
5331 open (iliptranpar,file=liptranname,status='old',action='read')
5332 call getenv_loc('TUBEPAR',tubename)
5333 open (itube,file=tubename,status='old',action='read')
5334 call getenv_loc('IONPAR',ionname)
5335 open (iion,file=ionname,status='old',action='read')
5337 ! print *,"Processor",myrank," opened file ISIDEP"
5338 ! print *,"Processor",myrank," opened parameter files"
5340 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5341 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5342 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5343 ! Get parameter filenames and open the parameter files.
5344 call getenv_loc('BONDPAR',bondname)
5345 open (ibond,file=bondname,status='old')
5346 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5347 open (ibond_nucl,file=bondname_nucl,status='old')
5349 call getenv_loc('THETPAR',thetname)
5350 open (ithep,file=thetname,status='old')
5351 call getenv_loc('ROTPAR',rotname)
5352 open (irotam,file=rotname,status='old')
5353 call getenv_loc('TORPAR',torname)
5354 open (itorp,file=torname,status='old')
5355 call getenv_loc('TORDPAR',tordname)
5356 open (itordp,file=tordname,status='old')
5357 call getenv_loc('SCCORPAR',sccorname)
5358 open (isccor,file=sccorname,status='old')
5359 call getenv_loc('FOURIER',fouriername)
5360 open (ifourier,file=fouriername,status='old')
5361 call getenv_loc('ELEPAR',elename)
5362 open (ielep,file=elename,status='old')
5363 call getenv_loc('SIDEPAR',sidename)
5364 open (isidep,file=sidename,status='old')
5366 open (ithep_nucl,file=thetname_nucl,status='old')
5367 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5368 open (irotam_nucl,file=rotname_nucl,status='old')
5369 call getenv_loc('TORPAR_NUCL',torname_nucl)
5370 open (itorp_nucl,file=torname_nucl,status='old')
5371 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5372 open (itordp_nucl,file=tordname_nucl,status='old')
5373 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5374 open (isidep_nucl,file=sidename_nucl,status='old')
5376 call getenv_loc('LIPTRANPAR',liptranname)
5377 open (iliptranpar,file=liptranname,status='old')
5378 call getenv_loc('TUBEPAR',tubename)
5379 open (itube,file=tubename,status='old')
5380 call getenv_loc('IONPAR',ionname)
5381 open (iion,file=ionname,status='old')
5382 call getenv_loc('IONPAR_NUCL',ionnuclname)
5383 open (iionnucl,file=ionnuclname,status='old')
5385 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5387 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5388 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5389 ! Get parameter filenames and open the parameter files.
5390 call getenv_loc('BONDPAR',bondname)
5391 open (ibond,file=bondname,status='old',action='read')
5392 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5393 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5394 call getenv_loc('THETPAR',thetname)
5395 open (ithep,file=thetname,status='old',action='read')
5396 call getenv_loc('ROTPAR',rotname)
5397 open (irotam,file=rotname,status='old',action='read')
5398 call getenv_loc('TORPAR',torname)
5399 open (itorp,file=torname,status='old',action='read')
5400 call getenv_loc('TORDPAR',tordname)
5401 open (itordp,file=tordname,status='old',action='read')
5402 call getenv_loc('SCCORPAR',sccorname)
5403 open (isccor,file=sccorname,status='old',action='read')
5405 call getenv_loc('THETPARPDB',thetname_pdb)
5406 print *,"thetname_pdb ",thetname_pdb
5407 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5408 print *,ithep_pdb," opened"
5410 call getenv_loc('FOURIER',fouriername)
5411 open (ifourier,file=fouriername,status='old',readonly)
5412 call getenv_loc('ELEPAR',elename)
5413 open (ielep,file=elename,status='old',readonly)
5414 call getenv_loc('SIDEPAR',sidename)
5415 open (isidep,file=sidename,status='old',readonly)
5417 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5418 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5419 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5420 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5421 call getenv_loc('TORPAR_NUCL',torname_nucl)
5422 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5423 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5424 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5425 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5426 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5427 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5428 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5429 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5430 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5431 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5432 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5433 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5434 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5437 call getenv_loc('LIPTRANPAR',liptranname)
5438 open (iliptranpar,file=liptranname,status='old',action='read')
5439 call getenv_loc('TUBEPAR',tubename)
5440 open (itube,file=tubename,status='old',action='read')
5441 call getenv_loc('IONPAR',ionname)
5442 open (iion,file=ionname,status='old',action='read')
5443 call getenv_loc('IONPAR_NUCL',ionnuclname)
5444 open (iionnucl,file=ionnuclname,status='old',action='read')
5447 call getenv_loc('ROTPARPDB',rotname_pdb)
5448 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5451 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5452 #if defined(WINIFL) || defined(WINPGI)
5453 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5454 #elif (defined CRAY) || (defined AIX)
5455 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5457 open (iscpp_nucl,file=scpname_nucl,status='old')
5459 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5464 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5465 ! Use -DOLDSCP to use hard-coded constants instead.
5467 call getenv_loc('SCPPAR',scpname)
5468 #if defined(WINIFL) || defined(WINPGI)
5469 open (iscpp,file=scpname,status='old',readonly,shared)
5470 #elif (defined CRAY) || (defined AIX)
5471 open (iscpp,file=scpname,status='old',action='read')
5473 open (iscpp,file=scpname,status='old')
5475 open (iscpp,file=scpname,status='old',action='read')
5478 call getenv_loc('PATTERN',patname)
5479 #if defined(WINIFL) || defined(WINPGI)
5480 open (icbase,file=patname,status='old',readonly,shared)
5481 #elif (defined CRAY) || (defined AIX)
5482 open (icbase,file=patname,status='old',action='read')
5484 open (icbase,file=patname,status='old')
5486 open (icbase,file=patname,status='old',action='read')
5489 ! Open output file only for CG processes
5490 ! print *,"Processor",myrank," fg_rank",fg_rank
5491 if (fg_rank.eq.0) then
5493 if (nodes.eq.1) then
5496 npos = dlog10(dfloat(nodes-1))+1
5498 if (npos.lt.3) npos=3
5499 write (liczba,'(i1)') npos
5500 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5502 write (liczba,form) me
5503 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5504 liczba(:ilen(liczba))
5505 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5507 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5509 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5510 liczba(:ilen(liczba))//'.mol2'
5511 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5512 liczba(:ilen(liczba))//'.stat'
5514 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5515 //liczba(:ilen(liczba))//'.stat')
5516 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5519 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5520 liczba(:ilen(liczba))//'.const'
5525 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5526 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5527 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5528 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5529 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5531 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5533 rest2name=prefix(:ilen(prefix))//'.rst'
5535 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5538 #if defined(AIX) || defined(PGI)
5539 if (me.eq.king .or. .not. out1file) &
5540 open(iout,file=outname,status='unknown')
5542 if (fg_rank.gt.0) then
5543 write (liczba,'(i3.3)') myrank/nfgtasks
5544 write (ll,'(bz,i3.3)') fg_rank
5545 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5550 open(igeom,file=intname,status='unknown',position='append')
5551 open(ipdb,file=pdbname,status='unknown')
5552 open(imol2,file=mol2name,status='unknown')
5553 open(istat,file=statname,status='unknown',position='append')
5555 !1out open(iout,file=outname,status='unknown')
5558 if (me.eq.king .or. .not.out1file) &
5559 open(iout,file=outname,status='unknown')
5561 if (fg_rank.gt.0) then
5562 write (liczba,'(i3.3)') myrank/nfgtasks
5563 write (ll,'(bz,i3.3)') fg_rank
5564 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5569 open(igeom,file=intname,status='unknown',access='append')
5570 open(ipdb,file=pdbname,status='unknown')
5571 open(imol2,file=mol2name,status='unknown')
5572 open(istat,file=statname,status='unknown',access='append')
5574 !1out open(iout,file=outname,status='unknown')
5577 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5578 csa_seed=prefix(:lenpre)//'.CSA.seed'
5579 csa_history=prefix(:lenpre)//'.CSA.history'
5580 csa_bank=prefix(:lenpre)//'.CSA.bank'
5581 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5582 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5583 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5584 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5585 csa_int=prefix(:lenpre)//'.int'
5586 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5587 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5588 csa_in=prefix(:lenpre)//'.CSA.in'
5589 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5592 write (iout,'(80(1h-))')
5593 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5594 write (iout,'(80(1h-))')
5595 write (iout,*) "Input file : ",&
5596 pref_orig(:ilen(pref_orig))//'.inp'
5597 write (iout,*) "Output file : ",&
5598 outname(:ilen(outname))
5600 write (iout,*) "Sidechain potential file : ",&
5601 sidename(:ilen(sidename))
5603 write (iout,*) "SCp potential file : ",&
5604 scpname(:ilen(scpname))
5606 write (iout,*) "Electrostatic potential file : ",&
5607 elename(:ilen(elename))
5608 write (iout,*) "Cumulant coefficient file : ",&
5609 fouriername(:ilen(fouriername))
5610 write (iout,*) "Torsional parameter file : ",&
5611 torname(:ilen(torname))
5612 write (iout,*) "Double torsional parameter file : ",&
5613 tordname(:ilen(tordname))
5614 write (iout,*) "SCCOR parameter file : ",&
5615 sccorname(:ilen(sccorname))
5616 write (iout,*) "Bond & inertia constant file : ",&
5617 bondname(:ilen(bondname))
5618 write (iout,*) "Bending parameter file : ",&
5619 thetname(:ilen(thetname))
5620 write (iout,*) "Rotamer parameter file : ",&
5621 rotname(:ilen(rotname))
5624 write (iout,*) "Thetpdb parameter file : ",&
5625 thetname_pdb(:ilen(thetname_pdb))
5628 write (iout,*) "Threading database : ",&
5629 patname(:ilen(patname))
5631 write (iout,*)" DIRTMP : ",&
5633 write (iout,'(80(1h-))')
5636 end subroutine openunits
5637 !-----------------------------------------------------------------------------
5640 use geometry_data, only: nres,dc
5642 ! implicit real*8 (a-h,o-z)
5643 ! include 'DIMENSIONS'
5644 ! include 'COMMON.CHAIN'
5645 ! include 'COMMON.IOUNITS'
5646 ! include 'COMMON.MD'
5649 ! real(kind=8) :: var,ene
5651 open(irest2,file=rest2name,status='unknown')
5652 read(irest2,*) totT,EK,potE,totE,t_bath
5655 ! AL 4/17/17: Now reading d_t(0,:) too
5657 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5660 ! AL 4/17/17: Now reading d_c(0,:) too
5662 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5665 read (irest2,*) iset
5669 end subroutine readrst
5670 !-----------------------------------------------------------------------------
5671 subroutine read_fragments
5675 use control_data, only:out1file
5678 ! implicit real*8 (a-h,o-z)
5679 ! include 'DIMENSIONS'
5683 ! include 'COMMON.SETUP'
5684 ! include 'COMMON.CHAIN'
5685 ! include 'COMMON.IOUNITS'
5686 ! include 'COMMON.MD'
5687 ! include 'COMMON.CONTROL'
5690 ! real(kind=8) :: var,ene
5692 read(inp,*) nset,nfrag,npair,nfrag_back
5694 !el from module energy
5695 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5696 if(.not.allocated(wfrag_back)) then
5697 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5698 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5700 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5701 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5703 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5704 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5707 if(me.eq.king.or..not.out1file) &
5708 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5709 " nfrag_back",nfrag_back
5711 read(inp,*) mset(iset)
5713 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5715 if(me.eq.king.or..not.out1file) &
5716 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5717 ifrag(2,i,iset), qinfrag(i,iset)
5720 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5722 if(me.eq.king.or..not.out1file) &
5723 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5724 ipair(2,i,iset), qinpair(i,iset)
5727 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5728 wfrag_back(3,i,iset),&
5729 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5730 if(me.eq.king.or..not.out1file) &
5731 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5732 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5736 end subroutine read_fragments
5737 !-----------------------------------------------------------------------------
5739 !-----------------------------------------------------------------------------
5743 ! implicit real*8 (a-h,o-z)
5744 ! include 'DIMENSIONS'
5745 ! include 'COMMON.CSA'
5746 ! include 'COMMON.BANK'
5747 ! include 'COMMON.IOUNITS'
5749 ! integer :: ntf,ik,iw_pdb
5750 ! real(kind=8) :: var,ene
5752 open(icsa_in,file=csa_in,status="old",err=100)
5753 read(icsa_in,*) nconf
5754 read(icsa_in,*) jstart,jend
5755 read(icsa_in,*) nstmax
5756 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5757 read(icsa_in,*) nran0,nran1,irr
5758 read(icsa_in,*) nseed
5759 read(icsa_in,*) ntotal,cut1,cut2
5760 read(icsa_in,*) estop
5761 read(icsa_in,*) icmax,irestart
5762 read(icsa_in,*) ntbankm,dele,difcut
5763 read(icsa_in,*) iref,rmscut,pnccut
5764 read(icsa_in,*) ndiff
5771 end subroutine csa_read
5772 !-----------------------------------------------------------------------------
5773 subroutine initial_write
5776 ! implicit real*8 (a-h,o-z)
5777 ! include 'DIMENSIONS'
5778 ! include 'COMMON.CSA'
5779 ! include 'COMMON.BANK'
5780 ! include 'COMMON.IOUNITS'
5782 ! integer :: ntf,ik,iw_pdb
5783 ! real(kind=8) :: var,ene
5785 open(icsa_seed,file=csa_seed,status="unknown")
5786 write(icsa_seed,*) "seed"
5788 #if defined(AIX) || defined(PGI)
5789 open(icsa_history,file=csa_history,status="unknown",&
5792 open(icsa_history,file=csa_history,status="unknown",&
5795 write(icsa_history,*) nconf
5796 write(icsa_history,*) jstart,jend
5797 write(icsa_history,*) nstmax
5798 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5799 write(icsa_history,*) nran0,nran1,irr
5800 write(icsa_history,*) nseed
5801 write(icsa_history,*) ntotal,cut1,cut2
5802 write(icsa_history,*) estop
5803 write(icsa_history,*) icmax,irestart
5804 write(icsa_history,*) ntbankm,dele,difcut
5805 write(icsa_history,*) iref,rmscut,pnccut
5806 write(icsa_history,*) ndiff
5808 write(icsa_history,*)
5811 open(icsa_bank1,file=csa_bank1,status="unknown")
5812 write(icsa_bank1,*) 0
5816 end subroutine initial_write
5817 !-----------------------------------------------------------------------------
5818 subroutine restart_write
5821 ! implicit real*8 (a-h,o-z)
5822 ! include 'DIMENSIONS'
5823 ! include 'COMMON.IOUNITS'
5824 ! include 'COMMON.CSA'
5825 ! include 'COMMON.BANK'
5827 ! integer :: ntf,ik,iw_pdb
5828 ! real(kind=8) :: var,ene
5830 #if defined(AIX) || defined(PGI)
5831 open(icsa_history,file=csa_history,position="append")
5833 open(icsa_history,file=csa_history,access="append")
5835 write(icsa_history,*)
5836 write(icsa_history,*) "This is restart"
5837 write(icsa_history,*)
5838 write(icsa_history,*) nconf
5839 write(icsa_history,*) jstart,jend
5840 write(icsa_history,*) nstmax
5841 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5842 write(icsa_history,*) nran0,nran1,irr
5843 write(icsa_history,*) nseed
5844 write(icsa_history,*) ntotal,cut1,cut2
5845 write(icsa_history,*) estop
5846 write(icsa_history,*) icmax,irestart
5847 write(icsa_history,*) ntbankm,dele,difcut
5848 write(icsa_history,*) iref,rmscut,pnccut
5849 write(icsa_history,*) ndiff
5850 write(icsa_history,*)
5851 write(icsa_history,*) "irestart is: ", irestart
5853 write(icsa_history,*)
5857 end subroutine restart_write
5858 !-----------------------------------------------------------------------------
5860 !-----------------------------------------------------------------------------
5861 subroutine write_pdb(npdb,titelloc,ee)
5863 ! implicit real*8 (a-h,o-z)
5864 ! include 'DIMENSIONS'
5865 ! include 'COMMON.IOUNITS'
5866 character(len=50) :: titelloc1
5867 character*(*) :: titelloc
5868 character(len=3) :: zahl
5869 character(len=5) :: liczba5
5871 integer :: npdb !,ilen
5875 ! real(kind=8) :: var,ene
5879 if (npdb.lt.1000) then
5880 call numstr(npdb,zahl)
5881 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5883 if (npdb.lt.10000) then
5884 write(liczba5,'(i1,i4)') 0,npdb
5886 write(liczba5,'(i5)') npdb
5888 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5890 call pdbout(ee,titelloc1,ipdb)
5893 end subroutine write_pdb
5894 !-----------------------------------------------------------------------------
5896 !-----------------------------------------------------------------------------
5897 subroutine write_thread_summary
5898 ! Thread the sequence through a database of known structures
5899 use control_data, only: refstr
5901 use energy_data, only: n_ene_comp
5903 ! implicit real*8 (a-h,o-z)
5904 ! include 'DIMENSIONS'
5906 use MPI_data !include 'COMMON.INFO'
5909 ! include 'COMMON.CONTROL'
5910 ! include 'COMMON.CHAIN'
5911 ! include 'COMMON.DBASE'
5912 ! include 'COMMON.INTERACT'
5913 ! include 'COMMON.VAR'
5914 ! include 'COMMON.THREAD'
5915 ! include 'COMMON.FFIELD'
5916 ! include 'COMMON.SBRIDGE'
5917 ! include 'COMMON.HEADER'
5918 ! include 'COMMON.NAMES'
5919 ! include 'COMMON.IOUNITS'
5920 ! include 'COMMON.TIME1'
5922 integer,dimension(maxthread) :: ip
5923 real(kind=8),dimension(0:n_ene) :: energia
5925 integer :: i,j,ii,jj,ipj,ik,kk,ist
5926 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5928 write (iout,'(30x,a/)') &
5929 ' *********** Summary threading statistics ************'
5930 write (iout,'(a)') 'Initial energies:'
5931 write (iout,'(a4,2x,a12,14a14,3a8)') &
5932 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5933 'RMSnat','NatCONT','NNCONT','RMS'
5934 ! Energy sort patterns
5939 enet=ener(n_ene-1,ip(i))
5942 if (ener(n_ene-1,ip(j)).lt.enet) then
5944 enet=ener(n_ene-1,ip(j))
5956 ist=nres_base(2,ii)+ipatt(2,i)
5958 energia(i)=ener0(kk,i)
5960 etot=ener0(n_ene_comp+1,i)
5961 rmsnat=ener0(n_ene_comp+2,i)
5962 rms=ener0(n_ene_comp+3,i)
5963 frac=ener0(n_ene_comp+4,i)
5964 frac_nn=ener0(n_ene_comp+5,i)
5967 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5968 i,str_nam(ii),ist+1,&
5969 (energia(print_order(kk)),kk=1,nprint_ene),&
5970 etot,rmsnat,frac,frac_nn,rms
5972 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5973 i,str_nam(ii),ist+1,&
5974 (energia(print_order(kk)),kk=1,nprint_ene),etot
5977 write (iout,'(//a)') 'Final energies:'
5978 write (iout,'(a4,2x,a12,17a14,3a8)') &
5979 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5980 'RMSnat','NatCONT','NNCONT','RMS'
5984 ist=nres_base(2,ii)+ipatt(2,i)
5986 energia(kk)=ener(kk,ik)
5988 etot=ener(n_ene_comp+1,i)
5989 rmsnat=ener(n_ene_comp+2,i)
5990 rms=ener(n_ene_comp+3,i)
5991 frac=ener(n_ene_comp+4,i)
5992 frac_nn=ener(n_ene_comp+5,i)
5993 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5994 i,str_nam(ii),ist+1,&
5995 (energia(print_order(kk)),kk=1,nprint_ene),&
5996 etot,rmsnat,frac,frac_nn,rms
5998 write (iout,'(/a/)') 'IEXAM array:'
5999 write (iout,'(i5)') nexcl
6001 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6003 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6004 'Max. time for threading step ',max_time_for_thread,&
6005 'Average time for threading step: ',ave_time_for_thread
6007 end subroutine write_thread_summary
6008 !-----------------------------------------------------------------------------
6009 subroutine write_stat_thread(ithread,ipattern,ist)
6011 use energy_data, only: n_ene_comp
6013 ! implicit real*8 (a-h,o-z)
6014 ! include "DIMENSIONS"
6015 ! include "COMMON.CONTROL"
6016 ! include "COMMON.IOUNITS"
6017 ! include "COMMON.THREAD"
6018 ! include "COMMON.FFIELD"
6019 ! include "COMMON.DBASE"
6020 ! include "COMMON.NAMES"
6021 real(kind=8),dimension(0:n_ene) :: energia
6023 integer :: ithread,ipattern,ist,i
6024 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6026 #if defined(AIX) || defined(PGI)
6027 open(istat,file=statname,position='append')
6029 open(istat,file=statname,access='append')
6032 energia(i)=ener(i,ithread)
6034 etot=ener(n_ene_comp+1,ithread)
6035 rmsnat=ener(n_ene_comp+2,ithread)
6036 rms=ener(n_ene_comp+3,ithread)
6037 frac=ener(n_ene_comp+4,ithread)
6038 frac_nn=ener(n_ene_comp+5,ithread)
6039 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6040 ithread,str_nam(ipattern),ist+1,&
6041 (energia(print_order(i)),i=1,nprint_ene),&
6042 etot,rmsnat,frac,frac_nn,rms
6045 end subroutine write_stat_thread
6046 !-----------------------------------------------------------------------------
6047 subroutine readpdb_template(k)
6048 ! Read the PDB file for read_constr_homology with read2sigma
6049 ! and convert the peptide geometry into virtual-chain geometry.
6051 ! include 'DIMENSIONS'
6052 ! include 'COMMON.LOCAL'
6053 ! include 'COMMON.VAR'
6054 ! include 'COMMON.CHAIN'
6055 ! include 'COMMON.INTERACT'
6056 ! include 'COMMON.IOUNITS'
6057 ! include 'COMMON.GEO'
6058 ! include 'COMMON.NAMES'
6059 ! include 'COMMON.CONTROL'
6060 ! include 'COMMON.FRAG'
6061 ! include 'COMMON.SETUP'
6062 use compare_data, only:nhfrag,nbfrag
6063 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6065 logical lprn /.false./,fail
6066 real(kind=8), dimension (3):: e1,e2,e3
6067 real(kind=8) :: dcj,efree_temp
6071 real(kind=8), dimension (3,20) :: sccor
6073 integer, dimension (:), allocatable :: iterter
6074 if(.not.allocated(iterter))allocate(iterter(nres))
6081 write (2,*) "UNRES_PDB",unres_pdb
6089 read (ipdbin,'(a80)',end=10) card
6090 if (card(:3).eq.'END') then
6092 else if (card(:3).eq.'TER') then
6095 itype(ires_old-1,1)=ntyp1
6096 iterter(ires_old-1)=1
6097 itype(ires_old,1)=ntyp1
6100 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6103 dc(j,ires)=sccor(j,iii)
6106 call sccenter(ires,iii,sccor)
6109 ! Fish out the ATOM cards.
6110 if (index(card(1:4),'ATOM').gt.0) then
6111 read (card(12:16),*) atom
6112 ! write (iout,*) "! ",atom," !",ires
6113 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6114 read (card(23:26),*) ires
6115 read (card(18:20),'(a3)') res
6116 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6117 ! & " ires_old",ires_old
6118 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6119 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6120 if (ires-ishift+ishift1.ne.ires_old) then
6121 ! Calculate the CM of the preceding residue.
6125 dc(j,ires_old)=sccor(j,iii)
6128 call sccenter(ires_old,iii,sccor)
6132 ! Start new residue.
6133 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6136 else if (ibeg.eq.1) then
6137 ! write (iout,*) "BEG ires",ires
6139 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6143 ires=ires-ishift+ishift1
6145 ! write (iout,*) "ishift",ishift," ires",ires,
6146 ! & " ires_old",ires_old
6147 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6149 else if (ibeg.eq.2) then
6151 ishift=-ires_old+ires-1
6153 ! write (iout,*) "New chain started",ires,ishift
6156 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6157 ires=ires-ishift+ishift1
6160 if (res.eq.'ACE' .or. res.eq.'NHE') then
6163 itype(ires,1)=rescode(ires,res,0,1)
6166 ires=ires-ishift+ishift1
6168 ! write (iout,*) "ires_old",ires_old," ires",ires
6169 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6172 ! write (2,*) "ires",ires," res ",res," ity",ity
6173 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6174 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6175 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6176 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6178 write (iout,'(2i3,2x,a,3f8.3)') &
6179 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6183 sccor(j,iii)=c(j,ires)
6185 if (ishift.ne.0) then
6186 ires_ca=ires+ishift-ishift1
6190 ! write (*,*) card(23:27),ires,itype(ires)
6191 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6192 atom.ne.'N' .and. atom.ne.'C' .and.&
6193 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6194 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6195 ! write (iout,*) "sidechain ",atom
6197 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6201 10 if(me.eq.king.or..not.out1file) &
6202 write (iout,'(a,i5)') ' Nres: ',ires
6203 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6207 ! write (iout,*) i,itype(i),itype(i+1)
6208 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6209 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6210 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6211 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6212 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6214 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6215 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6222 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6226 dcj=(c(j,i-2)-c(j,i-3))/2.0
6227 if (dcj.eq.0) dcj=1.23591524223
6232 else !itype(i+1).eq.ntyp1
6234 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6235 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6242 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6246 dcj=(c(j,i+3)-c(j,i+2))/2.0
6247 if (dcj.eq.0) dcj=1.23591524223
6252 endif !itype(i+1).eq.ntyp1
6253 endif !itype.eq.ntyp1
6255 ! Calculate the CM of the last side chain.
6258 dc(j,ires)=sccor(j,iii)
6261 call sccenter(ires,iii,sccor)
6265 if (itype(nres,1).ne.10) then
6269 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6270 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6277 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6281 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6282 if (dcj.eq.0) dcj=1.23591524223
6283 c(j,nres)=c(j,nres-1)+dcj
6284 c(j,2*nres)=c(j,nres)
6295 c(j,2*nres)=c(j,nres)
6297 if (itype(1,1).eq.ntyp1) then
6301 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6302 call refsys(2,3,4,e1,e2,e3,fail)
6309 c(j,1)=c(j,2)-1.9d0*e2(j)
6313 dcj=(c(j,4)-c(j,3))/2.0
6319 ! Copy the coordinates to reference coordinates
6325 ! Calculate internal coordinates.
6326 if (out_template_coord) then
6327 write (iout,'(/a)') &
6328 "Cartesian coordinates of the reference structure"
6329 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6330 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6332 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6333 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6334 (c(j,ires+nres),j=1,3)
6337 ! Calculate internal coordinates.
6338 call int_from_cart(.true.,out_template_coord)
6339 call sc_loc_geom(.false.)
6341 thetaref(i)=theta(i)
6346 dc(j,i)=c(j,i+1)-c(j,i)
6347 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6352 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6353 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6355 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6356 ! & vbld_inv(i+nres)
6361 cref(j,i+nres,1)=c(j,i+nres)
6371 end subroutine readpdb_template
6372 !-----------------------------------------------------------------------------
6374 !-----------------------------------------------------------------------------
6375 end module io_config